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The  image  algebra,  an  algebraic  structure  developed  by  G.  X.  Ritter  et.  al. 
specifically  to  meet  the  needs  of  digital  image  processing,  will  provide  the  mathemat- 
ical setting  for  this  investigation.  The  results  presented  in  this  dissertation  evolved 
as  a  direct  result  of  questions  that  arose  during  the  development  of  this  image  al- 
gebra. A  network  of  processors  for  a  parallel  computer  architecture  is  modeled  as  a 
subset  of  R'^.  An  image  is  represented  as  a  function  defined  on  such  a  subset.  The 
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subsets  of  concern,  here,  will  be  rectangular  arrays  with  m  rows  and  n  columns.  A 
template  will  be  represented  by  a  function  from  points  in  a  finite  rectangular  array 
into  images.  Basic  operands  and  operations  of  the  image  algebra,  such  as  addition, 
multiplication  and  convolution,  are  defined.  Relationships  between  image  algebra, 
matrix  (linear)  algebra,  and  polynomial  algebra  are  described.  The  two  main  ques- 
tions addressed  in  this  dissertation  are  the  problems  of  template  decomposition  and 
template  inversion. 

Methods  are  established  for  a  two-dimensional  shift  invariant  convolution  op- 
erator defined  on  a  rectangular  array  to  be  decomposed  into  sums  and  products  of 
operators  of  smaller  size.  The  main  result  of  this  type  is  that  5  x  5  shift  invariant 
operators  can  be  written  as  the  sum  and  product  of  at  most  five  3x3  shift  invariant 
operators.  The  second  is  that  a  5  x  5  operator  satisfying  certain  symmetry  condi- 
tions can  be  written  as  the  sum  and  product  of  at  most  three  3x3  operators.  The 
decomposition  methods  presented  here  are  suitable  for  a  machine  with  a  pipeline 
architecture. 

An  additional  focus  of  this  research  will  center  on  the  problem  of  extending 
results  valid  for  shift  invariant  operators  to  non-shift  invariant  operators.  In  par- 
ticular, necessary  and  sufficient  conditions  are  given  for  a  variant  template  to  be 
separable.  A  5  x  5  variant  template  can  always  be  written  as  the  sum  and  product 
of  at  most  seven  3x3  templates. 

The  mean  filter  with  five  nonzero  weights  all  equal  to  one  and  whose  configu- 
ration (or  shape)  is  in  the  form  of  a  cross  is  shown  to  be  invertible  when  defined  on 
an  m  x  m  array  if  and  only  if  neither  5  nor  6  divide  m. 
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CHAPTER  1 
INTRODUCTION 


This  dissertation  is  the  result  of  work  done  in  conjunction  with  an  investiga- 
tion of  tlie  structure  of  the  image  algebra,  an  algebraic  structure  for  use  in  image 
processing.  The  need  for  a  unifying  mathematical  theory  for  image  processing  algo- 
rithms specification,  and  a  more  powerful  basis  for  an  algebraically-based,  high  level 
programming  language,  led  to  the  development  of  the  image  algebra.  It  is  an  algebra 
that  operates  on  images,  and  whose  operators  reflect  the  types  of  transformations 
commonly  used  in  digital  image  processing  [2].  It  is  capable  of  expressing  all  types 
of  transformations  and  algorithms  commonly  used  in  image  processing,  in  terms  of 
its  operators  [28,38]. 

G.  Matheron  [18]  and  J.  Serra  [32]  initiated  the  use  of  image  algebra  in  image 
processing.  It  was  known  then  as  mathematical  morphology.  Sternberg  [33,34] 
also  investigated  it  independently.  This  algebra  was  based  on  the  operations  of 
Minkowski  addition  and  subtraction  of  sets  in  R'^  [24].  Mathematical  morphological 
is  about  the  study  of  shapes  and  patterns.  The  operations  of  Minkowski  are  often 
referred  to  as  dilation  and  erosion  operations.  Theses  morphological  operations 
techniques  can  be  applied  to  many  image  processing  and  image  analysis  problems. 
However,  they  cannot,  with  the  exceptions  of  very  few  simple  cases,  be  used  as  a 
basis  for  a  general  purpose  algebraically  based  language  for  digital  image  processing. 
Their  weakness  resides  in  the  fact  that  they  cannot  easily  express  many  commonly 
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used  transformations  such  as  linear  image  to  image  transformations  and  algorithms 
[19].  In  response  to  these  Hmitations,  and  in  order  to  establish  an  algebraic  structure 
capable  of  expressing  most  of  the  common  image  processing  algorithms,  G.  X.  Ritter 
at  the  University  of  Florida  began  the  development  of  a  more  general  image  algebra 
capable  of  expressing  all  linear  and  morphological  transformations  [23,25,26,28]. 

This  image  algebra,  will  provide  the  mathematical  setting  for  most  of  the  re- 
sults presented  in  this  dissertation.  In  Chapter  2  basic  operands  and  operations 
are  defined,  and  relationships  between  image  algebra  linear  algebra  and  polynomial 
algebra  are  discussed.  Isomorphisms  imbedding  the  linear  algebra,  and  the  poly- 
nomial algebra  into  the  image  algebra  give  rise  to  the  observation  that  a  linear  (or 
generalized)  convolution  is  equivalent  to  matrix  multiplication,  and  the  observation 
that  operator  inversion  (or  deconvolution)  is  equivalent  to  matrix  inversion.  In  the 
case  that  an  operator  is  circulant,  operator  inversion  is  equivalent  to  polynomial 
inversion.  In  the  usual  implementation  of  the  isomorphism  between  templates  de- 
fined on  rectangular  arrays  and  matrices,  a  shift  invariant  operator  corresponds  to  a 
block  Toeplitz  matrix  with  Toeplitz  blocks,  and  a  circulant  operator  corresponds  to 
a  block  circulant  matrix  with  circulant  blocks.  The  image  algebra  allows  us  to  give  a 
precise  formulation  of  the  decomposition  methods  for  non-shift  invariant  operators. 
Recall  that  some  frequently  used  non-shift-invariant  operators  include  the  Fourier 
Transform  or  the  Gabor  Transform  [7,11]. 

The  use  of  parallel  processing  has  come  to  play  an  increasingly  important  role 
in  image  processing  during  the  past  few  years.  Efficient  implementation  of  linear 
convolution  is  an  important  aspect  of  any  hardware  environment.  Some  classical  ex- 
amples of  two-dimensional  operators  include  the  Gaussian  mean  filter,  the  Laplacian 
filter,  and  the  Discrete  Fourier  Transform. 
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Template  decomposition  plays  a  fundamental  role  in  image  processing  algo- 
rithm optimization,  since  it  provides  means  for  reducing  the  cost  in  computation 
time,  and  therefore  increases  the  comptutational  efficiency  of  image  processing  al- 
gorithms. This  goal  can  be  achieved  in  two  ways,  either  by  reducing  the  number 
of  arithmetic  computations  in  an  algorithm  or  by  restructuring  the  algorithms  to 
match  the  structure  of  a  special  image  processing  architecture.  One  of  the  goals  of 
this  work  is  to  present  methods  for  the  decomposition  of  a  two-dimensional  shift- 
invariant  convolution  operator  in  a  form  suitable  for  a  parallel  image  processing 
machine  with  hardware  similar  to  the  PIPE  developed  at  the  Center  for  Manufac- 
turing Engineering  of  the  National  Bureau  of  Standards.  As  mentioned  by  O'Leary 
[22],  a  machine  of  this  type  is  able  to  efficiently  apply  a  3  x  3  convolution  operator 
to  a  256  X  256  image,  but  has  a  more  difficult  time  applying  a  5  x  5  operator. 
Another  example  of  a  pipeline  computer  is  ERIM's  CytoComputer  [34]  which  can 
deal  only  with  templates  of  size  3  x  3  or  smaller  on  each  pipeline  stage.  Thus,  un- 
less Fourier  Transform  methods  are  to  be  used,  operators  of  size  larger  than  3x3 
must  be  decomposed  into  sums  and  products  of  3  x  3  operators  to  fit  the  hardware 
restriction.  O'Leary  showed  that  if  an  operator  (viewed  as  a  matrix)  has  rank  n, 
then  it  can  be  written  as  the  sum  of  at  most  n  rank  1  (i.e.  separable)  operators. 
An  immediate  consequence  of  the  rank  one  method  described  by  O'Leary  in  [22]  is 
that  any  arbitrary  5x5  shift-invariant  operator  can  be  written  in  the  form 

t  =  (ti  e  t2)  +  (ts  ®  t4)  +  (t5  ®  te)  +  (t7  ®  tg)  +  (tg  ®  tio), 
where  each  tj  is  a  3  x  3  shift-invariant  operator  and  ®  denotes  the  convolution 
operation. 

While  the  results  on  shift  invariant  operators  are  presented  in  Chapter  3  as 
theorems  concerning  the  decomposition  of  polynomials  of  two  variables  into  sums 
and  products  of  lower  degree  polynomials,  the  theorems  discussed  in  Chapter  4  will 
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all  be  phrased  in  terms  of  sums  and  ®  operations  on  templates.  In  Chapter  3,  it  is 
proved  that  an  arbitrary  5x5  shift-invariant  template  t  defined  on  Z  x  Z  can  be 
decomposed  into  the  form 

t  =  (ti©t2)  +  (t3et4)  +  t5, 
where  each  tj  is  a  3  x  3  shift-invariant  template.  Thus,  the  number  of  3  x  3 
templates  needed  to  decompose  t  is  actually  half  of  that  indicated  by  the  rank 
method.  Example  5  of  Chapter  5  shows  that  this  decomposition  is  the  best  possible, 
in  the  sense  that  it  may  not  always  be  possible  to  decompose  an  operator  into  the 
sum  and  product  of  fewer  than  five  3x3  templates. 

The  most  commonly  used  template  in  the  decompositions  is  the  3  x  3  square 
template,  better  known  as  the  Moore  neighborhood.  It  has  the  following  form: 


Since  most  templates  used  in  image  processing  are  symmetric  or  skew  symmet- 
ric with  respect  to  the  x  or  the  y  axis,  a  discussion  of  the  problem  of  decomposing 
convolution  operators  exhibiting  some  type  of  symmetry  is  included  in  Chapter  3 
and  Chapter  4. 

An  second  approach  to  the  decomposition  of  two  dimensional  operators  based 
on  the  lower-upper  triangular  (LU)  factorization  of  a  rectangular  matrix,  will  result 
into  the  finite  sum  of  separable  operators.  A  comparison  of  the  polynomial  method 
and  the  LU  factorization  methods  is  given  in  Chapter  3  Section  4. 
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Results  concerning  the  decomposition  of  arbitrary  operators,  (i.e.  those  which 
may  not  be  shift  invariant)  are  presented  in  Chapter  4.  Necessary  and  sufficient 
conditions  for  a  variant  operator  to  be  separable,  i.e.  spHttable  into  the  convolution 
of  a  horizontal  operator  and  a  vertical  operator,  are  established.  Since  not  all 
templates  are  separable,  we  will  show  that  every  variant  template  can  be  written  as 
the  sum  of  separable  templates. 

We  conclude  by  pointing  out  that  template  decomposition  is  not  only  necessary 
in  the  case  where  a  special  image  processing  device  cannot  handle  large  templates 
but  is  also  desirable  when  the  problem  of  efficient  decomposition  is  concerned  in 
real-time  applications. 

Inverse  problems  have  come  to  play  a  central  role  in  modern  applied  mathemat- 
ics. While  classical  Hnear  and  non-linear  inverse  problems  abound  in  mathematical 
physics  [31],  they  are  also  important  in  such  imaging  areas  as  tomography  [31], 
remote  sensing  [31],  and  restoration  [30].  While  Rosenfeld  and  Kak  [30]  showed 
how  the  Wiener  filter  and  least  square  methods  can  be  used  to  restore  an  image, 
they  also  explained  the  equivalence  between  these  techniques  and  the  problem  of 
inverting  a  block  circulant  matrix  with  circulant  blocks.  In  1973,  G.  E.  Trapp  [36] 
showed  how  the  Discrete  Fourier  Transform  could  be  used  to  diagonalize  and  invert 
a  matrix  that  was  either  circulant  or  block  circulant  with  circulant  blocks.  While  the 
algebraic  relationship  between  circulant  matrices  and  polj'nomials  was  completely 
formulated  by  J.  P.  Davis  [6],  it  was  P.  D.  Gader  and  G.  X  .Ritter  [9,27]  who  made 
the  connection  between  polynomials  and  circulant  templates.  In  particular,  it  was 
found  that  a  circulant  template  defined  on  a  rectangular  array  with  m  rows  and 
n  columns  is  invertible  if  and  only  if  its  corresponding  polynomial  p{x,y)  has  the 
property  that  p(cj;^,u;,^)  ^  0,  for  all  0  <  j  <  n,  and  0  <  ^'  <  m.  (The  symbol  tUn 
denotes  the  root  of  the  unity,  exp(27rz77?).) 
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An  application  of  this  method  is  that  the  usual  3x3  mean  filter  defined  on  a 
rectangular  array  with  m  rows  and  n  columns  is  invertible  if  and  only  if  the  number  3 
does  not  divide  either  m  or  n.  Thus  for  example,  the  mean  filter  defined  on  an  array 
with  512  rows  and  512  columns  will  be  invertible,  while  the  mean  filter  defined  on 
an  array  with  512  rows  and  240  columns  will  not  be  invertible.  Gader  [8]  asked  if  a 
similar  result  can  be  obtained  for  the  von  Neumann  mean  filter.  The  von  Neumann 
mean  filter  is  the  template  defined  by 


t   = 


<1> 


The  main  result  of  Chapter  6  is  to  prove  that  if  t  is  defined  on  a  square  array 
with  m  rows  and  m  columns,  then  t  is  invertible  if  and  only  if  neither  5  nor  6  divide 


m. 


A  variety  of  examples  are  presented  in  chapter  4  to  both  illustrate  the  utility 
of  the  methods  developed  and  to  demonstrate  the  fact  that  some  theorems  are  best 
possible. 


CHAPTER  2 
IMAGE  ALGEBRA:  AN  OVERVIEW 


2.1.    Introduction 


The  following  chapter  contains  a  brief  review  of  the  fundamental  concepts  and 
notation  of  the  image  algebra.  We  will  need  these  ideas  to  extend  results  valid  for 
shift-invariant  operators  to  more  general  operators.  For  the  purpose  of  this  disser- 
tation, only  basic  definitions  and  operations  are  given.  The  relationship  between 
image  algebra  and  linear  algebra  will  also  be  described  in  this  chapter.  For  a  full 
and  detailed  description  of  all  image  algebra  operands  and  operations,  see  [29]. 

The  image  algebra  is  a  heterogeneous  algebraic  structure  specially  designed  to 
meet  the  needs  of  digital  image  processing.  It  is  used  as  a  model  and  tool  for  the 
development  of  local,  parallel  algorithms.  Several  commonly  used  image  processing 
transformations,  such  as  general  convolutions,  Discrete  Fourier  Transform,  and  edge 
detection  techniques  can  be  expressed  in  terms  of  the  image  algebra. 

The  use  of  image  algebra  in  digital  image  processing  was  initiated  by  Serra  [32], 
Miller  [19],  and  Sternberg  [33].  It  was  Ritter  [24]  who  showed  that  their  algebras 
were  all  equivalent  and  based  on  the  morphological  operations  of  Minkowski. 

An  image  algebra  is  an  algebra  whose  operands  are  images  and  subimages 
(or  neighborhoods).    It  deals  with  six  basic  type  of  operands,  namely,  value  sets, 
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coordinate  sets  (or  point  sets),  the  elements  of  the  value  and  the  coordinate  sets, 
images  and  templates. 

2.2.    Definitions  and  Background 

2.2.1.  Operands  of  the  Image  Algebra 

A  value  set  is  the  set  of  values  that  an  image  can  take.  It  can  be  any  semi-group 
with  a  zero  element.  The  most  commonly  used  ones  in  image  processing  are  the  set 
of  positive  integers  Z"*",  integers  Z,  rational  numbers  Q,  real  numbers  R,  positive 
real  numbers  R"*",  or  complex  numbers  C.  An  unspecified  value  set  will  be  denoted 
byF. 

A  coordinate  set  (or  a  point  set)  is  a  subset  of  an  n-dimensional  Euclidean  space, 
R*^,  for  some  n.  It  is  commonly  denoted  by  X  or  Y,  and  the  elements  of  such  sets 
are  denoted  by  lower  case  letters.  The  familiar  point  sets  are  the  rectangular  and 
hexagonal  arrays. 

Definition  2.1.   Let  X,  and  F  be  a  point  and  a  value  set,  respectively.  An  F  valued 
image  a  on  X  is  a  function  a  :  X  -^  F. 

Thus,  the  graph  of  an  F  valued  image  a  on  X  is  of  the  form 

a  =  {(x,a(x))  :  a(x)  G  F,  for  all  x  G  X}. 

The  set  X  is  called  the  set  of  image  points  of  a,  and  the  range  of  the  function  a 
is  called  the  set  of  image  values  of  a.  The  pair  (x,a(x))  is  called  a  picture  element 
or  a  pixel  ,  x  the  pixel  location  ,  and  a(x)  the  pixel  (or  gray)  value  .  We  will  denote 
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the  set  of  all  F  valued  images  on  X  by  F^^.    We  make  no  distinction  between  an 
image  and  its  graph. 

Definition  2.2.   An  image  a  :  X  — ^  F  has  finite  support  on  X  «/a(x)  ^  0  for  only 
a  finite  number  of  elements  x  G  X. 

Another  basic,  but  very  powerful  tool  of  the  image  algebra,  is  the  generalized 
template. 

Definition  2.3.    Let  X  and  Y  be  two  coordinate  sets,  and  let  F  be  a  value  set.  A 
generalized  F  valued  template  t  from  Y  to  X  is  a  function  t  :  Y  — »■  F     . 

Thus,  for  each  y  G  Y,t(y)  G  F^,  or  equivalently,  t(y)  is  an  F  valued  image 
on  X.  For  notational  convenience,  we  define  ty  =  t(y).  Thus, 

ty  =  {(x,ty(x)):xGX}. 

The  sets  Y  and  X  are  called  the  target  domain  and  range  space  of  t,  respectively. 
The  point  y  is  called  the  target  (or  domain)  point  of  the  template  t,  and  the  values 

ty(x)  are  called  the  weights  of  the  template  t  at  y.  Note  that  the  set  of  all  F  valued 

X  Y 

templates  from  Y  to  X  can  be  denoted  by  (F     )     . 

If  t  is  a  template  from  Y  to  X,  then  the  set 

S(ty)  =  {xGX:ty(x)  ^0} 

is  called  the  support  of  ty. 

If  t  is  an  F  valued  template  from  X  to  X,  and  X  is  a  subset  of  R",  then 
t  is  called  translation  invariant  (or  shift-invariant)  if  and  only  if  for  each  triple 
X,  y ,  z  G  X,  with  x  +  z  and  y  +  z  G  X,  we  have  that 
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ty(x)  =  ty+z(x  +  z). 

Note  that  a  translation  invariant  template  must  be  an  element  of  (F^)  .  Invariant 
operators  on  Z  x  Z  are  commonly  expressed  in  terms  of  polynomials  of  two  vari- 
ables. For  brevity,  we  will  frequently  refer  to  shift-invariant  templates  as  invariant 
templates.  A  template  which  is  not  necessarily  translation  invariant  is  called  trans- 
lation variant  or,  simply,  a  variant  template.  Translation  invariant  templates  occur 
naturally  in  digital  image  processing. 

2.2.2.  Operations  of  the  Image  Algebra 

The  basic  operations  on  and  between  F  valued  images  are  naturally  derived 
from  the  algebraic  structure  of  the  value  set  F. 

Let  X  be  a  subset  of  R".  Suppose  a  €  R^  and  t  €  (R"^)     . 
Addition  on  images  is  defined  as  follows:  If  a,  b  G  F     ,  then 

a  +  b  =  {(x,  c(x))  :  c(x)  =  a(x)  +  b(x),    x  G  X}. 

Higher  level  operations  are  the  ones  that  involve  operations  between  templates 

and  images,  and  between  templates  only. 

X  ^     , 

The  addition  of  two  templates  is  defined  pointwise.  If  s  and  t  E  (R     )     ,  then 

we  have 

(s  +  t)y(x)  =  Sy(x)  +  ty(x). 


Definition  2.4.    The  generalized  convolution  of  an  image  a  with  a  template  t  is 
defined  by 

a®t  =  {(y,b(y)):b(y)=   J]  a(x)ty(x),    y  G  X}. 

xgX 
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Linear  convolution  plays  a  fundamental  role  in  image  processing.  It  is  involved 
in  as  many  important  examples  as  the  Discrete  Fourier  Transform,  the  Laplacian, 
the  mean  or  average  filter  and  the  Gaussian  mean  filter. 

Definition  2.5.  // s  and  t  are  templates  on  X,  then  we  define  the  generalized 
convolution  of  the  two  templates  as  the  template  r  =  s  ©  t  hy  defining  each  image 
function  ry  hy  the  rule 

ry  =  {(z,ry(z))  :  ry(z)  =   ^  ty(x)sx(z),   where  z  €  X}. 

xgX 

Note  that  r  can  be  viewed  as  a  generalization  of  the  usual  notion  of  the  com- 
position of  two  convolution  operators.  If  the  templates  r  and  s  are  translation 
invariant,  then  except  for  values  near  the  boundary,  the  previous  definition  agrees 
with  the  usual  definition  of  polynomial  product. 

Note  also  that  if  s  and  t  are  two  invariant  templates,  then  r  would  be  an 
invariant  template  too.  Defining  r  at  an  arbitrary  y  €  Y  is  sufficient  to  define  the 
template  everywhere. 

Many  other  image  operations  are  described  in  detail  in  Ritter  et  al  [29].  A 
precise  investigation  of  the  linear  operator  ®  can  also  be  found  in  Gader  [8],  and  an 
extensive  study  of  other  non-linear  template  operations  can  be  found  in  Davidson 
[5]  and  Li  [16]. 

2.3.    Image  Algebra  and  Linear  Algebra 

If  X  is  a  finite  rectangular  subset  of  the  plane  with  m  rows  and  n  columns, 
then  it  can  be  linearly  ordered  left  to  right  and  row  by  row.  Thus,  we  can  write  X  = 
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{xi,X2,...,X7nn}-      Let      L^      denote     the     set      of     templates      on      X.      (i.e. 

Let  (Mr7m,+i*)  denote  the  ring  of  mn  x  mn  matrices  with  entries  from  F 
under  matrix  addition  and  multiplication.  For  any  template  t,  we  define  a  matrix 
Ml  =  {rriij)  where  rrnj  =  tx-(xj).  For  the  sake  of  notational  convenience  we  will 
write  tji  for  txj(xj). 

Define  the  mapping  V^  :  L^  ~^  ^mn  by  ^(t)  =  M^. 

The  next  Theorem  was  proved  by  Ritter  and  Gader  [8].  It  shows  that  there  is 
an  embedding  of  linear  algebra  into  image  algebra. 

Theorem  2.6.  The  mapping  V'  is  an  isomorphism  from  the  ring  {L^,+,(i))  onto 
the  ring  [Mmn, +■,*)■  That  is,  if  s,t  G  L-^,  then 

1.)  V'ls  +  t)  =  V'(s)  +  V'(t)  or  Ms+t  =  Ms  +  M^, 
2.)  V'(s  e  t)  =  V^(s)V'(t)  or  Mg^  t  =  MsMt, 
3.)  ijj  *^  one-to-one  and  onto. 

This  theorem  clearly  states  that  template  inversion  or  deconvolution  is  equiv- 
alent to  matrix  inversion.  Actually,  a  more  powerful  implication  of  this  theorem  is 
that  any  tool  available  in  linear  algebra  can  be  directly  applicable  to  problems  in 
image  algebra. 

Definition  2.7.  Let  'X.  be  an  m  x  n  coordinate  set.  We  say  that  the  mapping  (j)  : 
X  — >  X  is  a  circulant  translation  if  and  only  if  (j)  is  of  the  form 
^(x  +  h)  =  (x  +  h){mod{m,  n)),  for  some  h  G  X. 


-13- 

X  X 

Definition  2.8.    We  say  that  t  G  (F     )       ^5  circulant  if  and  only  if  for  every  cir- 

culant  translation  cj),  the  equation  tx(y)  =  ^d,(:x.)i^iy))  holds. 

This  last  definition  shows  that  a  circulant  template  is  completely  determined 
if  it  is  defined  at  only  one  point.  Translation  invariant  and  circulant  templates  are 
used  in  the  implementation  of  the  convolution  and  therefore,  are  directly  related  to 
the  Discrete  Fourier  Transform. 

For  the  purposes  of  this  dissertation,  X  will  usually  be  a  subset  of  Z  x  Z, 
where  Z  denotes  the  set  of  integers.  All  the  templates  we  will  consider  will  have 
finite  support  contained  in  a  set  of  the  form 

^m,n  =  WJ)  ■  \i\  <  m  and  |;|  <  n}. 

Therefore,  a  (2m  +  1)  x  (2n  +  1)  template  on  Z  x  Z  will  have  support  contained 
in  Xm,n  and  will  be  displayed  as  a  matrix  of  the  form 


-m,—n 


^0,-n 


im.  —  n 


i-m,0 


<  <0,0  > 


hnfi 


-m,n 


^0,n 


tm,n 


If  the  template  t  is  defined  as  above,  then  it  is  understood  that  all  weights 
tij  =  0,  for  all  {i,j)  such  that  |?"|  >  m  or  \j\  >  n.  The  element  <  ig.O  >  will  denote 
the  weight  at  the  center  pixel  location  (i.e.  tx{x)  =  io,o)-  Note  that  the  indexing 
here  differs  by  a  shift  from  the  one  given  above. 


CHAPTER  3 
SHIFT  INVARIANT  OPERATORS 


The  use  of  parallel  processing  has  been  increasing  for  the  past  years.  Linear 
convolution  is  widely  used  in  image  processing.  It  consists  of  applying  an  operation 
between  a  template  and  a  given  input  image,  pixelwise,  to  yield  an  output  image. 

One  of  the  main  reasons  for  template  decomposition  is  that  some  current  image 
processors  can  only  directly  implement  operations  on  very  small  templates  .  Most 
transforms  are  not  able  to  be  implemented  directly  on  a  parallel  architecture.  Con- 
sider, for  example,  a  parallel  image  processing  machine  with  hardware  similar  to 
the  PIPE,  developed  at  the  Center  for  Manufacturing  Engineering  of  the  National 
Bureau  of  Standards.  This  machine  is  capable  of  evaluating  a  small  template  ,  like 
a  3  X  3, but  cannot  easily  evaluate  a  template  of  larger  size. 

The  problem  of  template  decomposition  is  to  try  to  find  a  sequence  of  tem- 
plates ti,t2,...,t„  smaller  in  size  than  the  given  template  t,  such  that  when  we 
apply  an  operation  on  template  t  and  image  a,  it  will  be  equivalent  to  sequentially 
applying  operations  between  ti,t2,  ...,t„  to  the  image.  For  example,  in  the  case  of 
the  neighborhood  array  processor,  it  will  consist  of  decomposing  the  template  into  a 
sequence  of  factors,  where  each  factor  is  directly  implementable  on  the  architecture. 

The  motivation  behind  all  the  work  on  template  decomposition  is  to  increase 
the  speed  and  reduce  the  convolution  computation  time.  If  a  convolution  operator 
is  large  in  size,  then  the  cost  in  computation  time  can  be  prohibitive.  This  compu- 
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tational  complexity  can  be  reduced  if  the  given  template  in  decomposed  into  sums 
and  convolutions  of  smaller  size  templates. 

The  following  discussion  will  clarify  the  use  of  template  decomposition.    Sup- 
pose a  template  t  has  the  following  decomposition: 


t  =  imi^i)  +  (®j=?sj)- 


The  convolution  of  an  image  a  together  with  the  template  t  will  give 


®  t  —  a    ®    (®i=i''z'  ~(~  ®;=i^i) 


=  [(...((a®ri)©r2)...)®r/j]    +    [(...((a®  si)  ®  S2)...)  0  s;t] 

=  b  ®  c, 

where  b  =  [(...((a  ®  n)  ®  r2)...)  ®  rj^],  and  c  =  [(...((a®  si)  ®  S2)...)  ®  sj^]. 

When  a  ®  t  is  computed,  the  image  b  is  computed  first,  the  result  stored,  then 
the  image  c  is  computed  and  the  sum  of  the  two  is  taken. 

3.1.    Decomposition  Methods  for  General  Operators 

Templates  come  in  different  weights,  size  or  shape,  depending  on  the  image 
processing  application.  A  special  classs  of  invariant  templates  are  the  two  dimen- 
sional rectangular  templates.  The  support  of  a  rectangular  template  is  a  rectangular 
array. 
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While  a  5  X  5  operator  is  usually  defined  as  a  matrix  of  the  form 


^-2,2        ^-1,2 


-2,1 


•1,1 


^02 
^01 


^12 
^11 


^22 

^21 


^-2,0         ^-1,0        <  ^00  >        ^10  <20 

i-2,-1    ^-1,-1      h-i      h-i    h-i 

t-2-2      ^-1,-2         ^0,-2         ^1,-2      ^2,-2. 

we  will  express  all  the  theorems,  corollaries,  and  propositions  in  terms  of  polynomi- 
als. To  avoid  negative  exponents,  we  begin  all  indexing  with  the  integer  zero.  This 
implies,  if  n  =  m.  =  4,  that  the  remainder  term  r[x,y)  is  shifted,  so  that  its  value 
at  the  center  pixel  location  will  be  represented  by  r22-  Thus,  in  the  decomposition 
methods  presented  in  Sections  2,  3,  and  4,  it  will  never  be  necessary  to  shift  the 
image  before  applying  the  operator  r{x,y). 


Theorem  3.1.   Ift{x,y)  is  a  polynomial  defined  by  the  formula 


m      n 


f=0  j=0 
theii  there  exist  five  polynomials  pi{x,y),  P2{x,y),  qi{x,y),  92(2^)?/)  and  r[x,y)  such 
that 

t{x,  y)  =  pi(x,  y)qi(x,  y)  +  p2{x,  y)q2{x,  y)  +  r{x,  y), 


where 


2      2 


m—2  n—2 


<-\,:i 


Pl{x,y)  =  ^Y^piijx'y^ ,      qi{x,y)  =  J^  ^qujX-y, 
i=0  j=0  i=0  j=0 


2      2 


m-2  n-2 


P2{X:y)  =  ^'^P2ijx'y^,       q2{x,y)  =  X]  XI  ^'^ij-^'y^' 
i=0  j=0  i=0  j=0 
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and 


m—1  n  —  l 


r{x,y)=  Y^  ^rijx'y^. 
i=l  j=l 

Proof: 

If  a{y)  -  J2?=0^0iy^^  *h^"  ^y  *^^  Fundamental  Theorem  of  Algebra,  there  are 
two  polynomials  ai(y),  and  02(2/)  with  deg(Gi(?/))  <  2  and  deg(a2(?/))  <  n  -  2, 
such  that  a(y)  =  ai{y)a2{y).  Similarly,  if  c{y)  =  Y,?=oimiy\  then  there  are  two 
polynomials  ci(y),  and  C2{y)  with  deg(ci(y))  <  2  and  deg(c2(?/))  <  n  -  2,  such  that 
c(y)  =  ciiy)c2{y). 
If  we  choose 

pi(a;,j/)  =  ai{y)  +  x^ci{y), 
qiix^y)  ^  a2{y)  +  x'^-^C2{y), 
then  the  polynomial 

u(x,y)  =  i(x,j/)-pi(.T,y)9l(x,y) 
m      n 

z=0  ji=0 
will  have  the  property  that  if  Ujj  is  nonzero,  then  1  <  z  <  m  —  1. 

If    b{x)  =  X]^7^  UiQx\  then  the  polynomial  b{x)  can  be  decomposed  as  xbi{x), 
where  deg(6i(x)   <   m  -  2.     Similarly,  if  c?(.t)   =   Ei^l^  ^m^^  then  c?(x)  can  be 
decomposed  as  xdi{x),  where  deg{di{x))  <  m  —  2. 
Let 

P2{x,y)  =  x  +  xy'^,    and 
q2{^.y)  =  hi^)  +  di{x)y''-^. 


If 


r(x,y)  =  u{x,y)-p2{x,y)q2{x,y) 


m      n 


?=o  ?:=o 
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then  when  r^j  is  nonzero,  i  and  j  will  be  such  that  1  <  i  <  m  —  1,  and  1  <  j  <  n  —  3. 

Q.E.D. 

Corollary  3.2.   Ift{x,y)  is  a  polynomial  defined  by  the  formula 

4      4 

then  there  exist  five  polynomials  pi{x,y),  p2{x,y),  qi{x,y),  q2{x,y),  and  r{x,y)  such 
that 

t{x^  y)  =  Pi{x,  y)qi{x,  y)  +  P2{x,  y)q2ix,  y)  +  r(x,  y), 


where 


2      2  2       2 

pi{x,y)  =  YlY^P^ij^^y''^  qi(x^y)  =  ^^<iiij^'y-'^ 

i=0  j-0  i=0 j-0 

2       2  2       2 

i=0 j=0  i-0 j=0 


and 


3      3 

r{x,y)  =  ^^rijx'y^. 

Remark  3.3.  An  immediate  consequence  of  Theorem  3.1  is  that  if  t{x,y)  is  a  poly- 
nomial defined  by  t{x,y)  =  X^"— q  S?=o  ^ij^^y^ ■>  then  there  exist  7i  —  1  polynomials 
il{x,  y),  hi^,  2/)v,  tn-2{x,  y),  and  r{x,  y)  such  that 

t{x,y)  =  ti{x,y)  +  ...  +  tn-2{3.-,y)  +  r{x,y), 

where  each  ti{x,  y)  is  a  polynomial  which  can  be  written  as  a  product  of  polynomials 

of  the  form 

2      2 

P{x,y)  =  ^^Pijx'y^, 

i=0 j=0 
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and 

3      3 

i=lj-l 
3.2.    A  Characterization  of  Operators  Decomposable  into  a  Special  Form 

Theorem  3.4.   Let  t{x,y)  be  the  polynomial  defined  by  : 

m      n 

i=o j=0 

where  too,  ^Om  ^mO  and  tmn  are  all  nonzero. 

Let 

n  m 

n  rn 

(^{y) '^Y^miy\       and  d{x)  =  ^tiQx\ 

i=0  2  =  0 

Let  ai,Q'2,...,a„   be  the  roots  of  a{y);      (Si,  132,  ■■■,  l^m  be  the  roots  of  b{x); 
l\il2i  ■■■■•In  be  the  roots  of  c{y);      and  81,82,  ■■■,^m  be  the  roots  of  d{x). 
If  there  exists  an  ordering  of  the  roots  so  that  we  have  both  equations 

0102^1^2  =  I^llhl\l2i  and 
a3...an8s...8m  -  /?3...^m73---7n 
(or  both  equations 

aia2/3i(32  =  7172^1^2,  and 
a3...Q'n;53.../9m  =  73---7n^3---^m), 
satisfied,  then  there  exist  three  polynomials  p{x,y),  q{x,y)  and  r{x,y)  such  that 

t{x,  y)  =  p{x,  y)q{x,  y)  +  r{x,  y), 


where 
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2      2  m-2n-2 

p{x,y)  =  ^"^pijx'y^,       q{x,y)  =   J]  Z]  ^U^'?/^ 
i=Oj=0  i=0  j=0 


and 


m—l n— 1 
iz=l   j  =  i 


Conversely,  if 

t{x,y)  =  p{x,y)q{x,y)  +  r{x,y), 

where  p{x,y),q{x,y)  and  r{x,y)  are  defined  as  above,  then  there  exists  an  ordering 
of  the  roots  so  that  we  have  both  equations 

aia2^1^2  =  A^27l72i  «"^ 
a3...Qn^3---^m  =  /^3---A7?73---7n 
(or  both  equations 

Oi\OL2^\l^2  =  7172^1^2^  «"^ 
Q'3---ttn/^3---/^m  =  73---7n^3---<5m) 
satisfied. 
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Proof: 

Factor  a{y),b{x),c{y)  and  d[x)  to  obtain 

a{y)  =  toniy''  +  -  +  {t0l/t0n)y  +  (^OoAon)) 

=  tOniy'^  -  ("1  +  «2)2/  +  oi«2)(y""^  +  -  +  (-l)"a3...Qn),  (1) 

b{x)  =  tmnix'"  +  ■••  +  {tin/tmn)x  +  {ton/tmn)) 

=  (x2  -  (^1  +  p2)x  +  I3l/32){tmnx"'-'^  +  ...  +  (-Ij^^imn/^a-M,  (2) 
c{y)  =  imn(j/"  +  ...  +  {tml/tmn)y  +  (^mo/^mn)) 

=  {y^  -  (71  +  72)y  +  7l72)(^mn2/""^  +  •.•  +  (-l)"<m77,73-7n),  (3) 

c^(a;)  =  tmoix"^  +  ...  +  {tl0/tm0)x  +  (iooAmo)) 

=  <mO(a:^  -  (<5l  +  <52).T  +  <5i62)(x'"-2  +  ...  +  (-l)-63...^m),  (4) 

then  we  have  the  relations  : 

(-l)"Qia2...an  = /oo/'On,  (!') 

(-l)'"/?1^2-/?m  =  ton/tmn,  (2') 

(-l)"7l72---7n  =  tmO/tmn,    a.nd  (3') 

(-l)"^<5i,52...^r.j  =  <00Ar7.0.  (4') 

If  we  let 

o-liy)  =  Pll^2y'^  -  /5l/52(oi  +  a2)y  +  /?i/?2aifi2, 
and 

«2(y)  =  (-l)'"imn;53.../3^J/"-2  +  ...  +  (-l)'"+"imn/?3-Ana3...a„, 

then  by  substituting  formula  (2')  in  equation  (1),  it  will  be  true  that 

a{y)  =  a\{y)a2{y). 
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If 

di{x)  =  71723^^  -  7172(^1  +  h)^  +  7l72^1^2> 
and 

d2iy)  =  (-l)"i,„„73-7nx"^-2  +  ...  +  (-l)'"+"imn73-<^7,^3-^m, 

then  by  substituting  formula  (3')  into  equation  (4),  it  will  be  true  that 

d{x)  =  di{x)d2{x). 

Similarly,  we  can  rewrite  formulas  (2)  and  (3)  as 

b{x)  =  bi{x)b2{x)   and  c{y)  =  ci{y)c2{y). 

To  define  the  polynomial 

let  p(l,l)  =  0,  poi  be  the  coefficient  of  y^  in  ai(j/),  P2i  be  the  coefficient  of  y'  in 
ci{y),  piQ  be  the  coefficient  of  x^  in  di{x),  and  pj2  be  the  coefficient  of  x*  in  bi{x). 
To  define  the  polynomial 

let  qoj  be  the  coefficient  of  y^  in  02(^)5  Qm-2,j  be  the  coefficient  of  y^  in  C2(j/),  g/o 

be  the  coefficient  of  x^  in  c?2(a;),  and  qi^n-2  be  the  coefficient  of  x'^  in  62(^)5  ^'^^  set 

all  other  coefficients  equal  to  zero. 

If  we  set 

r{x,y)  =  t{x,y)-p{x,y)q{x,y) 
m      n 

i=0 j=0 
then,  it  follows  from  the  conditions  on  the  roots  that  if  rjj  is  nonzero,  then  1  <  ?  < 
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m  —  1  and  1  <  j  <  «  —  1. 

The  converse  can  be  proved  by  observing  that  if 

t{x,y)  =  p{x,y)q{x,y)  +  r{x,y), 
where 

^(^,y)  =  L/=i    Lj=i    ^ijx'yJ, 

then  equations  (l)-(4)  can  again  be  used  to  show  that  either 

aiO(2hh  =  A/^27172,  and 
a3...anSs...6m  =  (33...l3m^3---'fn 
or 

Qia2/5l/?2  =  7172^1^2  and 
a3---Otn/33--Pm  -  73. ..70^3... ^m- 


Q.E.D. 


Corollary  3.5.   Ift{x,y)  is  a  polynomial  defined  by  the  formula 

4      4 
i=0 j=0 

andt{x,y)  satifies  the  condition  of  Theorem  3.4,  then  there  exist  three  polynomials 
p{x,y),  q{x,y),  and  r{x,y)  such  that 

t{x,y)  =  p{x,y)q{x,y)  +  r{x,y), 


where 

2      2  2      2 

i=0 j=0  i=Q j=0 
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and 


3      3 


3.3.    Decomposition  Methods  for  Symmetric  Operators 

Most  operators  used  in  image  processing  are  symmetric  or  skew-symmetric 
with  respect  to  the  x  ov  y  axis.  The  following  section  deals  primarily  with  operators 
exhibiting  some  kind  of  symmetry. 

Definition  3.6.  A  polynomial  p{x)  =  anx^  +  a„_ix"~^  +  ...  +  aix  +  qq,  is  said 
to  satisfy  the  symmetric  property  with  respect  to  n,  if  gj  =  a^-i  for  all  i  =  0,  ...,n. 

Remark  3.7.  We  do  not  necessarily  assume  that  the  degree  of  p{x)  is  equal  to  n. 
For  example,  \i  p{x)  =  x,  then  p(x)  satisfies  the  symmetric  property  with  respect  to 
the  integer  2  . 

Definition  3.8.  A  polynomial  p{x)  is  said  to  satisfy  the  skew  symmetric  property 
with  respect  to  n,  if  Ui  =  —a^-i  for  all  i  =  0,  ...,n. 

Remark  3.9.  \{  p{x)  satisfies  the  symmetric  property  and  has  even  degree,  then  the 
exponents  of  x\  and  x'^"^  are  of  the  same  parity.  Therefore,  if  ,r  =  0  is  a  root  of 
multiplicity  A;,  then  p{x)  =  x^pi{x)^  where  p\{x)  is  of  even  degree. 
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ProPOSITION  3.10.    If  p{x)  is  a  polynomial  satisfying  the  symmetric,  or  the  skew 
symmetric  property  with  respect  to  some  integer  n,  then  a  non-zero  number  a  is  a 
root  ofp{x)  if  and  only  ifl/a  is  a  root  ofp{x). 

Proof: 

Since  ^(1/0-)  =  (±l/a")p(a)  and  q  is  a  non  zero  root,  the  polynomial  p{a)  = 
0  if  and  only  if  the  polynomial  p{llo)  =  0. 

Q.E.D. 

Proposition  3.11.  Ifp{x)  is  a  polynomial  of  even  degree  n  satisfying  the  symmet- 
ric property,  then  there  exist  two  numbers  ki  and  k2  such  that: 

p{x)  =  ax^^qi{x)q2{x)...qkj,x), 

where  a  is  the  leading  coefficient  ofp{x),  and  each  qj{x)  is  of  the  form  x    +  bjx  +  1. 

Proof: 

If  a  is  the  coefficient  of  the  highest  degree  of  x  in  p{x),  and  ki  is  the  largest 
integer  such  that  x^'-  divides  p{x),  then  p{x)  =  ax^'ipi(x),  where  pi{x)  satisfies  the 
symmetric  property.  Since  pi{x)  is  symmetric,  and  of  even  degree,  we  know  by 
Proposition  3.10  that  each  root  of  pi{x)  can  be  paired  with  its  reciprocal.  Thus, 
pi{x)  has  n/2  factors  of  the  form  [x  +  q){x  +  1/a)  =  x"^  -\- {a -\-  l/a)x  +  1.  Thus, 
pi{x)  =  (a-^  +  bix  +  1)...(.t2  +  6(„_^.3)/2-T  +  1)- 

Q.E.D. 

Proposition  3.12.  If  p{x)  =  Z1F=0"?'^*  '•^  ^  polynomial  satisfying  the  skew  sym- 
metric property  with  respect  to  an  even  integer  n,  then  there  exists  a  polynomial 
Pl{x)  such  that  p{x)  =  (x^  —  \)pi{x),  where  pi{x)  satisfies  the  symmetric  property. 
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PrOOF: 

Since  flj  =  — a„_j,  where  n  is  even, 

Similarly,  p(-l)  =  0.  Thus,  x^  -  1  is  a  factor  of  p{x). 

We  must  now  show  that  if  p(x)  =  (x^  -  l)pi{x),  then  pi{x)  satisfies  the  sym- 
metric property.  If 

pi{x)  =  6n-2a-""^  +  -  +  ^l--^-  +  ^0, 
then  note  that 


UQ 


=  -^0,^1  =  -6i,a„_i  =  6„_3,  and  «„  =  in-2- 


Since  p{x)  satisfies  the  skew  symmetric  property,  6„_2  =  ^0)  and  6n_3  =  6i. 
Note  also  that  aj,  =  bi^_2  —  hf,,  for  all  k  =  2,  ...,n  —  2. 

We  must  now  show  that  h).  =  6„_2_^^.  If  we  assume  inductively  that  6ju_2  = 
bf^_l,,  then,  since  p(a;)  satisfies  the  skew  symmetric  property, 

bk-2  -  bk  =  ak  =  -^n-k  =  -{K-k-2  -  K-k)-, 
so  that  6;;.  =  bn_2-k-  Thus,  pi{x)  is  symmetric. 

Q.E.D. 

Corollary  3.13.   If  p{x)  is  a  polynomial  satisfying  the  skew  symmetric  property 
with  respect  to  an  even  integer  n,  then  there  exist  two  numbers  ki  and  k2  such  that 

p{x)  =  a{x^  -  l)x^'qi{x)q2{x)...qk^{x), 

where  a  is  the  leading  coefficient  ofp{x),  and  qj{x)  is  of  the  form  (.t^  +  bjx  +  1),  for 
j  =  l,...,k2. 


-27- 
PrOOF: 

If  zero  is  a  root  of  p{x),  then  let  ki  be  its  multiplicity.  By  Proposition  3.12  and 
Remark  3.9,  there  is  a  symmetric  polynomial  pi{x),  such  that 

p{x)  =  anx'^'ix^  -  l)pi{x). 

Since  n  is  even,  and  pi(x)  satisfies  the  symmetric  property  with  respect  to  n, 
Proposition  3.11  implies  that  pi{x)  factors  into  a  product  of  quadratic  polynomials 
of  the  form  x    -\-  bx  -{-  I. 

Q.E.D. 

Remark  3.14.   If 

m      n 

i=0  j=0 
then  p{x,  y)  can  be  written  in  the  form 

The  polynomial  p{x,y)  can  also  be  written  in  the  form 


Definition  3.15.  Ifp{x,y)  =ViLoE]^otijx'y^  =  YJi=Qx'Pi{y)  -  then  the  poly- 
nomial  p{x,y)  is  said  to  satisfy  the  symmetric  (respectively  the  skew  symmetric) 
property  with  respect  to  y,  if  Pi{y)  satisfies  the  symmetric  (respectively  the  skew 
symmetric)  property,  for  all  i  =  0,  ...,n. 

A  similar  definition  can  be  given  for  the  variable  x. 
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Proposition  3.16.   Ifp{x,y)  is  any  polynomial,  then 

p{x,y)  =  Pi{x,y)  +  p2ix,y), 

where  pi{x,y)  satisfies  the  symmetric  property  with  respect  to  y,  andp2{x,y)  satisfies 
the  skew  sym.metric  property  with  respect  to  y. 

Proof: 

Let  pi(y)  =  Ei=o  "-"^^y^  and  p2{y)  =  Z]=o  "-"^^v' ■ 

Q.E.D. 

A  similar  proposition  can  be  given  for  tlie  variable  x. 

Proposition  3.17.  //  t{x,y)  =  ^^oY^]=otij^^y-'  has  the  property  that 
^00  =  ^On  =  ^mO  =  ^mn  =  0;  and  tiQ  =  Hqi,  tm,n-l  =  ktm-l,n>  ^m-1,0  =  ^'^ml) 
and  tQfi-i  —  ktifi,  where  k  =  +1  or  -1,  then  there  exist  three  polynomials  p{x^y), 
q{x,y),  and  r{x,y)  such  that  : 


t{x,  y)  =  p{x,  y)q{x,  y)  +  r(x,  y), 


where 


and 


2      2  77J-2  n-2 

p{x,  y)  =  ^^  Pijx'yK    <i{x,  y)=  J2Y1  ^ij^'y^' 

i=0  j=0  i=0  j=0 


771—1  n— 1 
H-T^,y)  =  Y^  ^rijx'yK 
i=\  j=l 
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Proof: 

If  k  =  1,  let 


p{x,  y)  =  X  +  y  +  xy'^  +  x'^y, 


q{x,y)  =  [toi+...  +  to^n-iy''~'^]  +  [hny'''^  +  -  +  im-l,nx'^    ^y"    "^1 

+  [tml^'^-^  +  -  +  im,n-2^"^"'j/""']  +  [^20-^  +  -  +  <m-2,0^"^-'], 
and 

r{x,  y)  =  t{x,  y)  -  p{x,  y)q{x,  y). 

li  k  =  -1,  let  p{x,y)  =  x  -  y  -  xy'^  +  .r^y,  q{x,y),  and  r(.T,y)  as  above. 

Q.E.D. 

Theorem  3.1  has  the  foolowing  corollaries 
Corollary  1  to  Theorem  3.1.  // 

m      n 

iz=0  j=0 
satisfies  the  symmetric  or  the  skew  symmetric  property  with  respect  to  y,  with  ^qo 

and  tmo  nonzero,  then  there  exist  three  polynomials  p{x,y),  q{x,y),  a7idr{x,y)  such 

that 

t{x,y)  =  p{x,y)q{x,y)  +  r{x,y), 


where 


am 


2       2  m-2n-2 

M^5 y)  =  ^^Pij^'y^j    9(^' y")^  JlYl ^ij^'y^^ 

i=0  j=0  ?:=o  j=0 

777—1  n— 1 

r{x,y)  =  5^  ^rijx'yK 

7=1    j  =  l 
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PrOOF: 

Note  that  since  t{x,y)  is  either  symmetric,  or  skew  symmetric  with  respect  to 
2/,  tQji,  and  tmn  are  also  both  nonzero.  Following  the  notation  used  in  the  proof  of 
Theorem  3.1,  the  polynomials  a(?/),  6(x),  c(y),  and  d{x)  are  such  that  b{x)  =  d{x), 
and  a{y)  and  c{y)  satisfy  the  symmetric  property.  Therefore,  by  Proposition  3.11 
or  Corollary  3.13,  each  root  of  a{y)  and  c{y)  can  be  paired  with  its  reciprocal,  and 
we  can  choose  ai{y)^a2{y),ci{y),  and  C2{y)  such  that  0102  =  1,  a^a^.-.an  =  1, 
7l72  =  1)  and  7374. ..7n  =  1.  Thus,  the  conditions  of  Theorem  3.1  are  satisfied,  and 
the  result  follows. 

Q.E.D. 
Corollary  2  to  Theorem  3.1.  // 

n       n 
i=0  j=0 

satisfies  the  symmetric  property  with  respect  to  both  x  and  y  and  if  the  matrix 
T  =  (tij)  is  symmetric,  then  there  exist  three  polynomials  p{x,y),  q{x,y)  andr{x,y) 
such  that  : 

t{x,  y)  =  p{x,  y)q{x,  y)  +  r{x,  y), 

where 

2       2  n-2n-2 

i=0  j=0  i=0  j=0 

and 

n—l n—l 

r{x,y)  =  Y^  ^rijx'y^. 

i=l  j=l 
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Moreover,  the  polynomials  p{x,y)^q{x^y),  and  r{x,y)  can  also  be  chosen  to  be  sym- 
metric with  respect  to  both  x  and  y. 

Proof: 

In  the  case  where  igo  ¥"  0'  the  polynomials  a{y),  b{x),  c{y)  and  d{x)  defined 
in  Theorem  3.1  all  have  the  same  coefficients  and  therefore  have  the  same  decom- 
positions. In  the  case,  where  ^oo  =  ^On  =  imO  =  imn  =  0,  the  result  follows  from 
Proposition  3.17. 

Q.E.D. 

3.4.     The  LU  Factorization 

For  a  detailed  discussion  of  how  the  LU  Factorization  Theorem  can  be  used  to 
implement  operators  in  such  a  way  that  their  operation  count  is  reduced,  see  Nikias 
et.  al  [21]. 

Proposition  3.18.  If  a  5x5  matrix  M  has  an  LU  factorization  (without  permu- 
tations), then  there  exist  two  rank  one  m.atrices  A  and  B  and  a  matrix  R  =  (?"jj) 
such  that 

M  =  A  +  B-{-  R, 

and    rij  =  r2j  =  rn  =  ri2  =  0,  for  ij  =  1,  ...5. 

Proof: 

Suppose  M  can  be  written  M  =  L  •  C/,  where  L  is  an  lower  triangular  matrix, 
and  U  is  an  upper  triangular  matrix.  Suppose  L  =  (/y),  and  U  —  (uij). 

Let  Ui  be  the  5x.5  zero  matrix  whose  r  row  is  the  i  row  of  U  and  let 
V  =  Us  +  U4  -\-  U^.  M  can  then  be  rewritten  as 


[1      ui2      Ui3      W14      ?«15] 
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M^L-U  =  L-{Ui  +  U2  +  V)  =  L-Ui  +  L-U2  +  L-V  =  A  +  B  +  R, 
where 

A  =  L-  Uu  B  =  L-U2,eind  R  =  L-V. 
The  matrix  A  is  a.  rank  one  matrix  of  the  form 

hi 
hi 
A=  hi 
hi 
hi. 

and  the  matrix  B  is  a  rank  one  matrix  of  the  form 

0 

h2 
B=      /32 

^52. 

On  the  other  hand,  i?  is  a  matrix  whose  first  two  rows  and  columns  are  all  zero. 

Q.E.D. 

Thus,  if  t  is  a  5  x  5  operator,  it  can  be  viewed  as  a  5  x  5  matrix.  If  t  has  an 
LU  factorization  without  permutations,  then  there  exist  two  rank  one  operators  u, 
and  V  and  an  operator  r  such  that  t  =  u  +  v  +  r. 

Since  rank  one  operators  are  separable  (i.e.  have  an  exact  factorization),  there 
exist  four  3x3  operators  ti,t2,t3,t4,  such  that 


[0      1      U23      U2A      "25. 


u  =  ti  ®  t2   and   v  =  13  0  t/[. 
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Therefore,  t  has  the  following  decomposition 

t  =  (ti®t2)  +  (t3®t4)  +  r, 

where 

<  ^11  >     tn     ^13 
^21  ^22     ^23 

^31  ^32      ^33 

While  Proposition  3.18  appears  to  be  as  good  as  the  results  in  Corollary  3.2, 
note  that  the  center  pixel  in  the  image  is  at  the  (1,1)  location.  Thus,  the  data  must 
be  shifted  before  it  can  be  convolved  with  the  operator  r. 

While  the  methods  discussed  in  Sections  1-4  give  techniques  for  the  decom- 
position of  an  n  x  n  operator,  most  operators  used  in  image  processing  are  either 
symmetric  or  skew  symmetric  with  respect  to  either  the  x  or  the  y  axis.  If  n  is  an 
odd  integer,  then  the  rank  of  the  operator  is  no  more  than  (n  -|- 1)/2  in  the  symmet- 
ric case,  and  no  more  than  (n  -  l)/2  in  the  skew  symmetric  case.  Using  the  rank 
method  or  the  LU  factorization  on  symmetric  or  skew  symmetric  operators  of  low 
rank,  it  is  frequently  the  case  that  a  decomposition  can  be  produced  which  will  be 
nearly  as  efficient  an  implementation  as  that  provided  by  our  methods.  Thus,  the 
methods  presented  here,  are  more  apt  to  provide  a  decomposition  superior  to  those 
provided  by  the  LU  factorization  on  templates  of  relatively  full  rank. 

3.5.    Grobner  Bases 

Grobner  bases  have  been  introduced  to  Bruno  Buchberger  [4]  in  his  disserta- 
tion in  1965.  They  are  named  after  W.  Grobner  [9],  his  thesis  advisor.  Buchberger's 


-34- 
method  of  Grobner  bases  is  a  technique  that  provides  algorithmic  solutions  to  a  va- 
riety of  problems  in  ideal  theory  as  well  as  solutions  to  numerous  problems  in  many 
other  fields.  It  is  a  general  purpose  method  for  multivariate  polynomial  computa- 
tions. It  is  implemented  in  all  major  computer  algebra  systems  such  as  MACSYMA, 
MAPLE,  MuMATH,  REDUCE,  SCRATCHPAD  II,  etc...  In  spite  of  its  simplicity, 
the  Buchberger  algorithm  solves  a  wide  range  of  problems  from  symbolic  algebra  to 
computational  geometry. 

In  this  dissertation,  we  will  apply  it  to  systems  of  algebraic  equations  and  linear 
homogeneous  equations  with  polynomial  coefficients. 

3.5.1.  Terminology 

Let  F  be  a  field,  and  let  F[yi,y2, --^yn]  be  the  ring  of  n-variate  polynomials  over 
F.  The  letters  f,g,  h  will  stand  as  typed  variables  for  polynomials  in  F[y\,y2,  ■■■,  yn]] 
H,  G  for  finite  subsets  of  F[j/i,  y-2,  ...,yn];t,u  for  power  products  of  the  form  y\^  ...J/n"; 
a,b  for  elements  in  the  field  F,  and  i,j  for  natural  numbers. 

Definition  3.19.    If  H  is  defined  as  above,  then  the  ideal  generated  by  H  is  the  set 

Ideal{H)  =  {Y^fihi  :  hi  G  HJi  G  F[yi,2/2, -,  2/n]} 
of  all  linear  combinations  of  the  elements  of  H  in  F[yi,  y25  ■•■•.2/rj]- 

Definition  3.20.  A  total  ordering  <  on  the  power  products  y-^  ...yn  is  called  ad- 
missible if  1  =  y\---yn  i^  minimal  under  <  and  multiplication  by  a  power  product 
preserves  the  ordering. 
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Examples  for  such  orderings  are  the  total  degree  ordering  (i.e.  1  <  yi  <  2/2  < 
y?  <  yi2/2  <  2/2  <  2/1  <  •••  ^^  ^^^  bivariate  case),  and  the  purely  lexicographic 
ordering  (i.e.  1  <  J/i  <  y?  <  •■■  <  y2  <  yiy2  <  vlvi  <  •••  <  yf  <  2/l?/2  <  ■••  ^^  ^^^ 
bivariate  case). 

Given  a  fixed  ordering  <,  we  will  denote  the  coefficient  of  t  in  g  by  coefF(5r,  i),  the 
leading  power  product  of  g  with  respect  to  <  by  lpp{g),  and  the  leading  coefficient 
of  g  with  respect  to  <  by  lc(^). 

Definition  3.21.  Given  a  nonzero  polynomial  f,  and  two  other  polynomials  g,  and 
h  in  F[yi,y2,...,yn],  we  say  that  g  reduces  to  h  modulo  f  or  g  -^h  f^  '/  ^^d  only 
if  there  exist  u,b,  b  ^  0,  such  that  b.u.lpp{h)  is  identical  with  some  monomial  of  g, 
and  f  —  g  —  b.u.h. 

Such  a  reduction  can  be  viewed  as  a  generalized  division,  deleting  a  monomial 
in  g.  The  result  /  is  always  strictly  smaller  than  g,  with  respect  to  the  ordering  <  . 

Example  1.  Let  F  =  Q,  g  =  y\y2y?>  -  y?,  and  h  =  yiys  -  y2.  If  we  let  u  =  y2,  and 
6=1,  and  we  use  the  lexicographic  ordering,  then  f  =  g  —  b.u.h  =  y^—  y-^. 

Definition  3.22.  We  say  that  f  is  in  normal  form,  (or  reduced  form)  if  and  only 
if  there  is  no  /'  such  that  f  -^h  f  ■ 

Example  2.  Let  g  and  h  be  as  in  Example  1.  If  /^i  =  6,  and  /i2  =  yiVl  -  Vl^  then 
/  =  yl  —  y^  of  Example  1  is  irreducible  modulo  {/?i,  112] ■ 

Another  reduction  of  ^  is  the  following  one:  g  ^  f'  =  yiy^-Ui  ^  .f"  =  y2-yi- 
In  this  case,  /i"  is  irreducible  modulo  {/?i,  /i2}- 
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Therefore,  both  /  and  /"  are  normal  forms  of  g  modulo  {/^i,  /?2}- 
Hence,  there  exist  several  reduction  paths  that  can  lead  to  a  normal  form, 
modulo  a  basis  (or  a  set)  H. 

Definition  3.23.    The  set  H  is  called  a  Grobner  basis  if  and  only  if  each  g  has  a 
unique  normal  form  modulo  H. 

We  will  note  that  different  orderings  will  give  rise  to  different  Grobner  bases, 
and  that  a  Grobner  basis  is  not  necessarily  unique.  It  is  possible  to  reduce  a  poly- 
nomial of  the  basis  modulo  another  polynomial  of  the  basis,  and  hence  obtain  a 
smaller  basis. 

As  a  matter  of  experience,  using  total  degree  ordering  yields  much  better  com- 
puting time  in  general.  On  the  other  hand,  when  using  the  purely  lexicographic 
ordering,  the  resulting  basis  is  in  triangular  form  (i.e.  each  polynomial  introduces 
at  most  one  new  variable). 

3.5.2.  Application  to  Systems  of  Equations  and  Algorithm 

Given  a  system  of  algebraic  equations  H^  the  Grobner  basis  can  be  used  to: 

-decide  whether  H  is  solvable; 

-decide  whether  H  has  finitely  or  infinitely  many  solutions; 

-find  all  solutions  of  H,  in  the  finite  case. 
The  following  algorithm  shows  how  the  Grobner  basis  method  using  the  lexico- 
graphic ordering  can  be  applied  to  determine  the  answer  to  a  system  of  equations. 
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Algorithm 

G  =  Grobner  basis  of  H. 

liieG 

then  H  is  unsolvable. 

Else  if  there  exist  an  integer  j  such  that  no  leading  power  product  of  the 
polynomials  in  G  is  of  the  form  y^,  for  some  e, 
then  H  has  infinitely  many  solutions, 
else  do  /=  the  (only)  polynomial  in  G  H  F[yi], 
Xi  =  {{a):f{a)  =  0}, 
for  i  =  l,...,n  —  1, 

^i+1  =  0, 

for  all  {ai,...,a{)  G  Xi,  do 

K  =  {g{ai,...,ai,yi+i)  :  g  G  G  H  F[yi,  ...,yi^i] 

\F[yi,...:yi]}, 
f  =  g.c.d.  of  the  polynomials  in  K, 
Xi+i  =  Xi  U  {(ai,...ai,a)  :  /(a)  =  0}. 
Finally,  Xn  contains  all  the  solutions  of  H. 

The  Grobner  basis  method  was  used  in  some  examples  of  Chapter  5  to  show 
that  a  given  template  cannot  be  decomposed  into  a  certain  form.  Several  examples 
are  presented  in  Chapter  5.  The  equations  were  set  up  as  a  system  which  contains 
both  linear  and  nonlinear  equations.  For  example,  for  a  5  x  5  template,  we  will 
have  25  equations,  17  of  which  are  not  linear.  To  avoid  having  to  use  Grobner  bases 
theory  to  solve  a  system  of  nonlinear  equations,  we  reduced  the  problem  to  roots 
of  polynomials  of  one  variable  to  decompose  the  boundary  of  the  template.  Once 
that  outer  boundary  has  been  matched  we  only  face  a  system  of  linear  equations  to 
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match  the  second  boundary.  While  the  Grobner  bases  method  can  be  very  slow  in 
solving  nonlinear  equations,  it  does  a  very  fine  job  on  solving  linear  equations.  All 
computations  were  done  on  MAPLE. 


CHAPTER  4 
VARIANT  OPERATORS 


We  now  turn  our  attention  from  results  concerning  shift-invariant  templates 
to  analogous  results  for  templates,  which  are  not  necessarily  shift-invariant.  Recall 
that  the  Discrete  Fourier  Transform  and  the  Gabor  Transform  are  examples  of  vari- 
ant templates.  We  will  try  to  determine  the  class  of  templates  that  are  separable 
(i.e.  decomposable  into  the  convolution  of  a  horizontal  template  and  a  vertical  tem- 
plate or  vice- versa).  Since  not  all  templates  are  separable,  we  will  show  that  every 
template  can  be  written  as  the  sum  of  separable  templates. 

4.1.    Decomposition  of  General  Operators 

Definition  3.1.  A  template  t  defined  on  ZxZ  has  column  rank  one  if  it  has  finite 
support  and  there  exists  a  column  vector  u  such  that  for  each  pixel  location  x  in 
ZxZ  and  every  column  vector  v  in  t{x),  there  exists  a  scalar  A  such  that  v  =  Au. 

A  similar  definition  could  be  given  for  row  rank  one  templates. 

Definition  4.2.  A  template  t  defined  on  Z  x  Z  and  having  finite  support  is  said 
to  be  horizontal  if,  for  every  x  in  ZxZ,  there  exists  a  unique  nonzero  row  vector  v 
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in  t(x).  Similarly,  a  template  t  is  said  to  be  vertical  if,  for  every  x  in  Z  x  Z,  there 
exists  a  unique  nonzero  column  vector  v  in  t(x). 

Theorem  4.3.  //t  is  a  column  rank  one  template  defined  on  Z  x  Z  and  has  finite 
support,  then  there  exist  an  invariant  vertical  template  t^  and  a  (possibly  variant) 
horizontal  template  t2  such  that  t  =  ti  ®  t2. 

Proof: 

For  every  x  in  Z  x  Z,  t(x)  is  an  m  x  n  array.  This  array  contains  n  column 
vectors  Vi,  V2, ...,  Vn-  Since  the  template  has  rank  one,  for  each  i  =  1,...,?^,  there 
exists  a  scalar  Aj  and  a  fixed  column  vector  u  such  that  Vj  =  A^u.  Therefore, 

t(x)  =  [u]  •  [Ai    A2    ...    An]      (outer  product), 

where  [u]  is  an  invariant  vertical  template  and  [Ai  A2  ...  A^]  is  a  (possibly  variant) 
horizontal  template. 

Q.E.D. 

Remark  4.4.  Let  X  be  an  77  x  ra  subset  of  Z  x  Z,  and  let  F  be  a  value  set.  The 
two-dimensional  Fourier  template  K  is  defined  by 

K{u,v)  =  —exp{  —  {uXi  +  vYi)2-Ki/n}, 
n 

where  Xi.Yi  G  F^  are  defined  by 

A'l  =  {{{x,y),z)  -.x^z}  and  Y'l  =  {{{x,y),z)  ■.y  =  z}. 

Even  though  it  can  be  shown  that  K  is  a  separable  variant  template,  (i.e.  it 
splits  into  the  product  of  a  horizontal  and  a  vertical),  K  has  neither  row  rank  one 
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nor  column  rank  one.  If  the  Fourier  template  K  did  have  column  rank  one,  then  by 
Theorem  4.3,  it  could  be  factored  into  the  product  of  an  invariant  vertical  template 
and  a  variant  horizontal  template. 

Definition  4.5.  If  a  template  t  defined  on  Z  X  Z  has  finite  support,  then  we  say 
that  t  has  column  rank  at  most  k  if  there  exist  k  column  vectors  ui,U2,  ...,uj^  such 
that  for  every  pixel  location  x  in  Z  x  Z  and  every  column  vector  v  in  t(x),  there 
exist  k  scalars  A^,  A2, ...,  A;;,  such  that  v  =  J2i=l  '^i^i- 

Q.E.D. 

Theorem  4.6.  If  a  column  rank  k  template  t  defined  on  Z  x  Z  has  finite  support, 
then  there  exist  k  column  rank  one  (possibly  variant)  templates  ti,  t2, ...,  tj^  such 
thatt  =  ti  +  t2  +  ...  +  t]^. 

Proof: 

By  Definition  4.5,  for  every  pixel  location  x  in  Z  x  Z,  and  every  column  vector 
Vj  in  t(x),  there  exist  k  scalars  Afi,  Aj2,  •■•,  ^ik  s^ch  that  Vf  =  J2j=l  ^ij^j- 

Let  ^rij  =  \ijUj  for  i  =  l,...,n.  If  tj(x)  =  [vij  V2j  ...  v^j]  for  j  =  l,...,k, 
then  each  tj  is  a  column  rank  one  template  and  t  =  ti  +  t2  +  ...  +  ty^.. 

Corollary  4.7.  //  a  template  t  defined  on  Z  X  Z  has  finite  support  and  column 
rank  k,  then  there  exist  k  vertical  invariant  templates  ri,  r2, ...,  r^^^  C'^^  ^  horizontal 
variant  templates  si,S2, .-.,  Sjt,  such  that 

k 
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ProoF: 

This  result  follows  immediately  from  Theorems  4.3  and  4.6. 

Q.E.D. 

Note  that  a  result  similar  to  Corollary  4.7  can  be  proved  for  row  rank  k  tem- 
plates. 

Theorem  4.8.  Let  t  he  a  template  on  a  finite  rectangular  subset  of  the  plane  with 
m  rows  and  n  columns,  and  let  M^  be  the  matrix  described  in  section  4  of  Chapter 
2.  If  Ars  denotes  the  matrix  [m„Y,._i)+2  „(j_i)4s]  fori  =  l,...,n  and  j  =  l,...,n, 
then  the  tem.plate  t  is  separable  if  and  only  if  each  matrix  Aj-s  has  rank  one  for  all 
r  =  1,  ...,m  and  s  —  1,  ...,m. 

Proof: 

Since  for  every  r  =  l,...,m  and  s  =  l,...,r77,  Ars  is  an  n  x  n  rank  one  matrix, 
there  exist  a  column  vector  v^s  and  a  row  vector  h^s  such  that  Mrs  =  v^shrs- 

Let  the  column  vectors  vis,  V2s, ...,  v^g  form  the  s  n  x  n  block  of  the  block 
diagonal  matrix  V,  for  5  =  1,  ...,m. 

On  the  other  hand,  let  the  i^/j  elements  of  the  n  row  vectors  h^i,  hr2?  •■•i  hm 
form  the  diagonal  elements  of  the  {r,i)  diagonal  block  of  the  block  matrix  H,  for 
r  =  1, ...,  m. 

A  direct  computation  V  H  will  give  the  matrix  M^.  Under  the  isomorphism  tp, 
V  and  H  represent  the  image  of  a  vertical  and  a  horizontal  template,  respectively. 
We  then  conclude  by  Theorem  2.6,  that  t  is  separable. 

Q.E.D. 
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COROLLARY  4.9.  Let  t  he  a  template  on  a  finite  rectangular  subset  of  the  plane  with 
m  rows  and  n  columns,  and  let  M^  be  the  matrix  described  in  the  previous  paragraph. 
IfArs  denotes  the  matrix  [mj,(^r-l)+iMJ-'^)+s'^  for  i  =  l,...,n  and  j  =  I,  ...,n,  then 
the  template  t  can  be  written  as  the  sum  of  k  separable  templates  if  and  only  if  each 
matrix  Ars  has  rank  at  most  k  for  all  r  =  1, ...,  m  and  5  =  1, ...,  m. 

Proof: 

Since  for  every  r  =  l,...,m  and  s  =  l,...,m,  Ars  is  an  n  x  n  matrix  of 
at  most  rank  k,  there  exist  k  column  vectors  Vjls,  v^^, ...,  v^'^  and  k  row  vectors 
hl^,  h2„  ...,  h^,  such  that  if  aIs  =  vishis.  for  i  =  1, ...,  ^-,  then  Ars  =  A}.,  +  /l^,  + 
...  +  Ars-  For  j  =  1, ...,  A',  Ais  is  an  n  x  n  rank  one  matrix,  and  therefore  Theorem 
4.9  applies.  If  for  every  j  =  1,...,A;  Mj  is  the  block  matrix  [Ais],  for  r  =  l,...,m 
and  5  =  1, ...,  m,  then  there  exists  two  matrices  Vj  and  Hj  such  that  A^g  =  Vj  ■  Hj. 

Since  Mt  =  Ml  +  M2  +  ...  +  Mk  impHes  that  M^  =  Vi- Hi  +  V2- H2  +  ...  +  Vk- Hk, 
Theorem  4.8  can  now  be  used  to  conclude  that  the  template  t  is  the  sum  of  k 
separable  templates. 

Q.E.D. 

Remark  4.10.  Similarly,  if  we  let  Ars  denote  the  matrix  [m„(-j_i)_,_,.  „(s_i)^.j],  then 
we  can  expect  the  same  type  of  results  as  in  Theorem  4.8.  The  only  difference  is  that 
the  horizontal  and  the  vertical  templates  are  convolved  in  the  reverse  order.  Recall 
also  that  the  two-dimensional  Fourier  template  is  separable.  It  is  straightforward  to 
check  that  the  Fourier  template  has  the  property  that  each  of  the  matrices  Ars  has 
rank  one. 

Remark  4.11.   Note  that  according  to  the  definition  of  the  convolution  of  two  tem- 
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plates  t  ©s,  it  is  as  if  we  are  "sliding"  the  first  template  t  along  the  second  one 
s.  Therefore,  in  order  to  minimize  the  number  of  equations  and  therefore  minimize 
the  complexity  of  the  solution,  it  will  be  very  "helpful"  to  have  the  first  template 
t  invariant.  The  goal  of  the  next  five  theorems  is  to  make  this  last  statement  more 
precise.  In  fact,  in  this  case,  polynomials  (rather  than  template  notation)  can  be 
used  in  the  decomposition. 

Theorem  4.12.    If  t  is  a  5  x  5  variant  template,  then  there  exist  seven  templates 
^15 ^2, t3, t4, t5, tg,  and  tj  such  that 

t  =  (ti  ®  t2)  +  (t3  ®  t4)  +  (t5  ®  tg)  +  t7, 

where  ti,t3,  and  ts  are  3x3  invariant  templates  and  t2,t4,t6,  and  tj  are  3x3 
variant  templates. 


Proof: 

Let  t  be  of  the  form 


t  = 


^-2,-2  ^-2,-1  ^-2,0  ^-2,1  ^-2,2 

^-1,-2  ^-1,-1  ^-1,0  ^-1,1  ^-1,2 

^0,-2  ^0,-1  <  ^0,0  >  ^0,1  <0,2 

^1,-2  ^1,-1  ^1,0  ^1,1  h,2 

^2,-2  ^2,-1  ^2,0  ^2,1  ^2,2 


and  let  us  define  the  following  shift  invariant  templates: 


■  1  = 


1 

0 

0" 

0 

<0> 

0 

t3  = 

0 

0 

1. 
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0        0         1 

0  <0  >     0 

1  0         0 


,     and  t5 


0  1         0 

1  <  0>     1 
0         1         0 


In  addition,  let  us  define  the  following  variant  templates: 


t2         = 


t-2-2      ^-2,-1        0 
i_l_2      <0>      <i,2 
0  <2,1         ^2,2 


t4 


0  t_2,l       ^-2,2 

il  _2      <  0  >      i_i,2 
^2,-2       ^2,-1  0 


and 


0  i_2,0  0 

^0 -2      <  0  >      ^0,2 

0  ^2,0  0 


Note  that  a  routine  calculation  using  the  definition  of  the  ®  operation  can  be 
used  to  show  that  the  template  ty  =  t  -  (ti  ®  t2  +  ta  ®  t4  +  ts  ®  te)  is  a  3  x  3 
template.  Thus,  t  can  be  decomposed  as  the  sum  and  product  of  at  most  seven 
3x3  templates. 

Q.E.D. 
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4.2.    Decomposition  of  Symmetric  Operators 

Theorem  4.13.  If  t  is  a  (2m  +  1)  x  (2n  +  1)  variant  template  symmetric  with 
respect  to  both  axes,  then  there  exist  four  templates  ti,t2,t3,  and  t4,  such  that 

t  =  (ti®t2)+t3  +  t4, 

where  tj  is  a  3  x  3  invariant  template,  t2  is  an  (2m  - 1)  x  (2n  -  1)  variant  templates 
symmetric  with  respect  to  both  axes,  13  is  a  (2m  +  1)  x  1  vertical  template,  and  14 
is  a  1  X  (2n  +  1)  horizontal  template. 

Proof: 

Let  t2  be  a  {2rn  —  1)  x  (2n  —  1)  template  symmetric  with  respect  to  both  axes, 
whose  left  top  m.  x  n  corner  is  given  by  the  equations: 


t2      ^-      ^    = 

^  — m+z,  — n+j 


t_,„+,-_„+,-     if     (z,j)  =  (0,0),(0,l),(l,0),(l,l) 


r(«,i) 


otherwise, 


where 


r(^j)  =  { 


^-m-\-i-n-\-j  +  ^-m-\-i-n-^j-2  "  t  <  J 

^-m^i-n-\-j  +  ^-m-\-i-2-n-\-j  ^^  ^  >  J 
^-m-{-i,-n+j  +  ^-rn+i-2,—n+j 

.   +t-m+i-n+j-2  +  ^-m+i-2-n+j-2  if  ^  =  J- 


Because  of  the  symmetry  of  t2,  it  is  easy  to  know  all  the  other  entries  of  that 
template. 


Now  if  ti  is  the  template  tj 
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1          0  -1 

0       <  0  >  0 

-1         0  1 


,  then  t  -  (ti  et2)  is 


simply  the  sum  of  a  vertical  template  t^  and  a  horizontal  template  14. 


Q.E.D. 


Theorem  4.14.   If  t  is  a  (2m  +  1)  x  (2n  +  1)   variant  template  symmetric  with 
respect  to  both  axes,  then  there  exist  five  templates  ti,  13,13, 14,  and  t^  such  that 

t  =  (ti®l2)  +  (l3®t4)+l5, 

where  ti  and  t^  are  3x3  invariant  templates  and  l2,l4?  c^^d  Is  are  (2m  —  1)  x 
(2n  —  1)  variant  templates  symmetric  with  respect  to  both  axes. 

Proof: 

The  main  idea  in  this  proof  is  to  decompose  the  cruciform  template  s  (i.e. 
S"  =  0,  whenever  i  ^  0  ov  j  ^  0)  derived  from  the  proof  of  Theorem  4.14. 


If 


t3  = 


0  1         0 

1  <  0>     1 
0         1         0 


and  I4  is  the  cruciform  template  of  size  (2m  -  1)  x  (2??  -  1)  with  same  extremal 
elements  as  in  s,  then  a  simple  computation  shows  that  1  —  (li  ®  I2  +  I3  ®  ^4)  is  a 
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(2m  +  1)  X  (2n  +  1)  template  symmetric  with  respect  to  both  axes. 


Q.E.D. 


Theorem  4.15.  7/t  is  a  (2m  +  l)  x  (2n  +  l)  variant  template  skew-symmetric  with 
respect  to  both  the  vertical  and  the  horizontal  axes,  then  there  exist  two  templates  ti 
and  t2  such  that 

t  =  ti  ©t2, 

where  ti  is  the  invariant  template  defined  by 


ti  = 


1  0  -1 

0       <0  >       0 

-1         0  1 


and    t2  is  a  (2m  —  1)  x  (2n  —  1)  variant  template  symmetric  with  respect  to  both 


axes. 


Proof: 

If  t2  =  [t2ij]  is  the  (2m  —  1)  x  {'2n  —  1)  template  symmetric  with  respect  to 
both  axes  defined  in  the  proof  of  Theorem  4.13,  then  it  is  a  matter  of  computation 
to  check  that  ti  ®  t2  is  the  given  (2m  +  1)  x  (2??  +  1)  skew-symmetric  template  t. 

Q.E.D. 

Theorem  4.16.   If  t  is  a  {2m  +  1)  x  (2m  +  1)  totally  symmetric  template  and  if 
the  corner  values  are  zero,  then  there  exist  three  templates  ti,t2,  and  t^  such  that 


t  =  (tiet2)  +  t3, 
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where  ti  is  a  3  x  3  invariant  template  and  t2  and  t-^  are  (2m  -  1)  x  {2m  -  1)  totally 
symmetric  variant  templates. 

Proof: 

If  the  corner  values  of  t  are  zero,  then  each  boundary  has  either  (2m  -  1)  or 
(2n  -  1)  elements  which  form  the  boundaries  of  a  (2m  -  1)  x  {2n  -  1)  template  t2 
symmetric  with  respect  to  both  axes. 


If  ti 


0  1         0 

1  <  0  >     1 
0         1         0 


then  the  template  ta  =  t  -  ti  ®  t2  is  a  {2m  -  l)x 


(2n  -  1)  template  symmetric  with  respect  to  both  axes. 


Q.E.D. 


CHAPTER  5 


EXAMPLES 


In  the  examples  presented  below,  we  have  indicated  the  value  in  the  center 
pixel  location  by  <  x  >  .  Except  for  Example  6,  all  templates  are  shift-invariant. 

Example  1. 
The  template 


-13  2  7  2  -13 

2  17  22  17       2 

7  22  <  27  >  22       7 

2  17  22  17       2 

-13  2  7  2  -13 


(used  by  Haralick  [4])  satifies  the  conditions  of  Corollary  2  to  Theorem  3.1.  Using 
the  techniques  of  Corollary  2  it  has  the  following  decomposition 


-13  -19.73  -13 

-19.73     <  -28.9433  >     -19.73 

-13  -19.73  -13 


© 


1  -1.672  1 

-1.672     <  1.5411  >     -1.672 
1  -1.672  1 


+  fS.39' 
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Using  the  rank  method  (note  that  h  has  rank  2),  the  template  h  can  be  de- 
composed as 


h 


+ 


1  1.5275  1   ■ 

1.5282  <  2.3191  >  1.5282 

1  1.5275  1   . 

1  1.3333  1 

1.2396  <  1.6528  >  1.2396 

1  1.3333  1 


e 


-13  19.8575  -13 

21.7364     <  -33.2024  >     21.7364 

-13  19.8575  -13 


® 


<  14.5208  > 
2 


Using  the  LU  factorization,  the  template  h  can  be  decomposed  into 


h  = 


1      1.5182      1 

1.5182  <  2.3049  >  1.5182 

1      1.5182      1 


® 


-13      21.7364      -13 
21.7364  <  -36.3433  >  21.7364 
L  -13      21.7364      -13 


+ 


17.3077  23.0769  17.3077 
23.0769  <  30.7692  >  23.0769 
17.3077         23.0769  17.3077 


Thus,  the  first  decomposition  provides  a  more  efficient  implementation  of  h 
than  the  other  two  decompositions. 
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Example  2. 

Using  Corollary  3.2,  the  template 


-5  -4045 
-8  -10  0  10  8 
10  -20  <  0  >  20  10 
-8  -10  0  10  8 
-5-4  0  4       5 


can  be  decomposed  into  the  form 

-5  -4  -5 

s=      -8     <  -10  >     -8 

-5  -4  -5 


1         0         -1 

0  <0  >       0 

1  0         -1 


+  [-12     <0>     12]. 


Note  that  the  template  s  is  skew  symmetric  with  respect  to  the  y-axis  and 
symmetric  with  respect  to  the  x-axis. 


Example  3. 

The  Laplace  template 
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0  0  -1  0  0 

0  0         0  0  0 

-1  0  <4>  0  -1 

0  0         0  0  0 

0  0  -1  0  0 


satifies  the  conditions  of  Proposition  3.6.  It  can  be  decomposed  into 


1 

0 

Q- 

0 

<  0> 

0 

e 

0 

0 

1. 

0  0-1 

0      <0>      0 
-10  0 


+  [4; 
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Example  4. 
The  template 


0     1         0  10 

0     0         0  0     0 

0     0  <  0  >  0     0 

10         0  0     1 

0     0         0  0     0 


cannot  be  decomposed  into  the  form  (ti®t2)  +  ts,  where  each  tj  is  a  3  x  3  operator. 
Using  the  LU  factorization  or  the  rank  method,  t  can  still  be  decomposed  as  the 
sum  and  product  of  four  3x3  operators.  Note  that  this  template  is  symmetric  with 
respect  to  the  y-axis.  Thus,  the  hypothesis  in  Corollary  1  to  Theorem  3.1  that  the 
corner  values  are  different  from  zero  is  necessary. 
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Example  5. 
The  template 


t  = 


0     1         0  10 

0     0         2  0     0 

0     1  <  1  >  1     0 

11         0  11 

0     0         0  0     0 


cannot  be  decomposed  into  the  form  (ti  ®  t2)  +  (ta®  14),  where  each  tj  is  a  shift- 
invariant  3x3  operator.  Note  that  this  template  is  symmetric  with  respect  to  the 
y-axis. 
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Example  6. 

Let  X  =  {(1,1),  (1,2),  (2,1),  (2,2)},  and  let  t  be  the  variant  template  defined 
by  the  rules 


t(l,l)  = 


1     1 
0     0 


t(l,2)  = 


0  0 

1  0 


t(2,i: 


1 

0' 

,andt(2,2)  = 

"0 

0 

0 

0 

1 

1 

The  matrix  M^  associated  with  t  is  given  by 


Mf 


110  0 
0  0  10 
10  0  0 
0     0     11 


Note  that  while  the  template  t  has  rank  one  for  all  x  in  X,  it  fails  to  be 
separable. 


CHAPTER  6 
OPERATOR  INVERSION 


6.1.    Introduction 


Inverse  problems  have  come  to  play  a  central  role  in  modern  applied  mathemat- 
ics, such  as  in  mathematical  physics,  in  imaging  areas  such  as  tomography,  remote 
sensing,  and  restoration.  Rosenfeld  and  Kak  [30]  explained  the  equivalence  between 
techniques  such  as  the  Wiener  filter  and  the  least  square  methods  and  the  problem 
of  inverting  a  block  circulant  matrix  with  circulant  blocks.  In  1973,  G.E.Trapp  [35] 
showed  how  the  Discrete  Fourier  Transform  could  be  used  to  diagonalize  and  invert 
a  matrix  that  was  either  circulant  or  block  circulant  with  circulant  blocks.  While  the 
algebraic  relationship  between  circulant  matrices  and  polynomials  was  completely 
formulated  by  J. P.  Davis  [6],  it  was  P.D.Gader  and  G.X.Ritter  [9,27]  who  made 
the  connection  between  polynomials  and  circulant  templates.  A  circulant  template 
defined  on  a  rectangular  array  with  m  rows  and  n  columns  is  invertible  if  and  only 
if  its  corresponding  polynomial  p{x,y)  has  the  property  that  p(a;^,a;f^)  7^  0,  for 
a-11  0  <  j  <  n,  and  0  <  A:  <  m.  (The  symbol  a>„  denotes  the  root  of  the  unity, 
exp(27rz/n).) 

An  application  of  this  method  is  that  the  usual  3x3  mean  filter  defined  on  a 
rectangular  array  with  m  rows  and  n  columns  is  invertible  if  and  only  if  the  number 
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3  does  not  divide  either  m  or  n.   Gader  asked  if  a  similar  result  applies  to  the  von 
Neumann  mean  filter.  The  von  Neumann  mean  filter  is  the  template  t  defined  by 


t   = 


<1> 


The  main  result  of  this  chapter  is  to  prove  that  if  t  is  defined  on  a  square  array 
with  m  rows  and  m  columns,  then  t  is  invertible  if  and  only  if  neither  5  nor  6  divide 
m. 

In  this  section,  we  will  briefly  describe  the  relationship  between  polynomial 
algebra  and  image  algebra.  More  precisely,  the  relationship  between  circulant  tem- 
plates and  quotient  rings  of  polynomial  rings.  Important  results,  such  as  how  to 
decompose  circulant  templates  and  how  to  obtain  minimal  local  decompositions  of 
separable  circulant  templates  can  be  found  in  Gader  [8]. 

6.2.    The  von  Neumann  Template 

6.2.1.  Preliminaries 


Definition  6.1.   A  m  x  m  matrix  C  —  c{i,j)  is  said  to  be  circulant  if  and  only  if 
for  every  5  e  Z,    Cij  =  c^i+s)(rnod   m),(j+s){mod  m)- 
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Thus,  C  is  a  matrix  of  the  form 


CO 

ci 

C2   ■ 

•  c,„_  1 

C,71-l 

CO 

Ci   . 

•   c,„_2 

C2 

C3 

C4   . 

ci 

c\ 

C2 

C3   • 

CO 

Each  row  is  just  the  previous  row  cycled  forward  one  step,  so  that  the  entries 
in  each  row  are  just  a  cyclic  permutation  of  those  in  the  first  row.  We  will  write 
C  =  circ(co,ci,...,c,„_i). 


The  m  x  m  permutation  matrix  P  = 


0 

1 

0   . 

.  0 

0 

0 

1    . 

.  0 

0 

0 

0   . 

..    1 

1 

0 

0  . 

..  0 

is  called  the  basic 


permutation  matrix. 

A  circulant  matrix  C  =  circ(co,ci,  ...,c^,_i),  can  be  associated  with  the  poly- 
nomial pc{x)  =  CO  +  cix  +  ...  +  Cm-ix"^'^.  If  C  =  circ(co,ci,...,c,„_i),  then 
C  -  pci.P)-  Because  of  this  representation,  circulant  matrices  have  a  very  nice 
structure  which  can  be  related  to  P.  Since  a  matrix  is  circulant  if  and  only  if 
it  commutes  with  the  matrix  P,  the  product  of  two  circulants  is  again  a  circulant. 
Thus,  the  set  Cm  of  all  m  x  m  circulant  matrices  is  a  subring  of  the  set  of  all  matrices 
of  order  m.  Furthermore,  circulant  matrices  commute  under  multiplication. 

Let  R[x]  denote  the  ring  of  polynomials  in  one  variable  with  coefficients  in  R, 
and  R[.r]/(a:'"  -  1)  the  quotient  ring  of  polynomials  modulo  x'"  -  1.  Multiplication 
of  elements  of  this  ring  are  performed  by  replacing  x'"  by  1. 
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If  we  define  the  mapping  ??  :  C,„  -^  R[x]/{x"^  -  1),  by  r]{C)  =  Pc{x),  then  7/ 
is  a  ring  isomorphism  [6]. 

Definition  6.2.  Let  B  he  an  mn  x  mn  matrix.  We  say  that  B  is  block  circulant 
with  circulant  blocks  of  type  {m,n)  if  and  only  if  there  exist  nxn  circulant  matrices 
Co,Ci,...,Cm-l  such  that  B  -  circ(Co,Ci,  ...,Cm-l) 

As  in  the  circulant  case,  given  a  block  circulant  matrix  of  the  form 
C  =  circ(co,ci,...,c,„n_i),  with  circulant  blocks  of  type  (m,7z),  we  can  associate 
a  polynomial  pc{x,y)  in  two  variables  as  follows: 

Pc{x,y)  =  (co  +  ciy  +  ...  +  c„_iy"-^)  +  x(c„  +  c^+iy  +  ...  +  csn-iy""^)  + 
...  +  x"^-l(c(^_i)„  +  c^ni-l)n+\y  +  -  +  Cmn-lJ/"~^)- 


As  in  the  one  variable  case,  R[x,y]/{x^  -  l,y"  -  1)  denotes  the  quotient  ring 
of  polynomials  modulo  x'"  —  1,  and  y"  —  1. 

If  we  define  the  mapping  ^  :  Cmn  ^  Ti[x,y]/{x'^-l,y''-l)hy  e{C)  =  pc{x,y), 
then  0  is  a  ring  isomorphism  [6]. 

Therefore,  given  a  circulant  matrix  M  e  Cmn,  M  is  invertible  if  and  only  if 
^(M)  =  PMi^^y)  is  a  unit  in  Ii[x,y]/{x'^  -  l,y"  -  1). 

Let  Um  =  exp(-27rij/m),  where  i  =  \/^,  and  let  A  be  the  diagonal  matrix 
diag(u;4),  for  J  =  0,l,...,77z  -  1.  It  can  be  proved  that  the  permutation  matrix 
described  above,  satisfies  the  equation  P  =  FAF*,  where  F  is  the  one  dimensional 
Discrete  Fourier  Transform  of  order  m,  and  *  denotes  the  conjugate  transpose.  Since 
F  is  a  unitary  matrix,  we  have  that  F*  =  F^^ .  Therefore,  if  C  is  circulant,  then 

C  =  pc{P)  =  PciFAF*)  =  Fpc{h)F*  =  FAF*, 
where  A  is  the  diagonal  matrix,  diag(pc('^m)),  for  J  =  0, 1, ...,  777  -  1,  whose  diagonal 
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terms  are  the  Fourier  coefficients.   Thus,  any  circulant  matrix  can  be  diagonahzed 
by  the  Fourier  matrix. 

Similarly,  if  F  is  the  mn  x  mn  two-dimensional  Discrete  Fourier  matrix,  and 
M  is  a  block  circulant  matrix  with  circulant  blocks,  then  we  can  write  M  =  PDF*, 
where  F  is  the  mn  x  mn  two-dimensional  Fourier  matrix,  and  D  is  a  diagonal  matrix. 

Therefore,  if  the  Fourier  coefficients  of  M  are  all  nonzero,  then  D  is  invertible 
and  the  inverse  of  M  is  given  by  M"^  =  FD~^F*  [6]. 

Let  t  be  an  invariant  circulant  template.  The  following  theorem  has  been 
proved  in  Gader  [8]  . 

Theorem  6.3.   Let  t  he  a  circulant  template  on  the  coordinate  set  X.    The  corre- 
sponding matrix  M^  is  a  block  circulant  matrix  with  circulant  blocks. 

Since  we  can  write  Mt  =  FDF*,  where  D  =  diag(pt(^m,'^n))  ^or 
j  =  0,...,m-l,  and  k  =  0, ...,  n-1,  M^  is  invertible  if  and  only  if  p^{u^rn,^n)  ¥"  0^ 
for  every  _;  and  k. 

Example:  Let  t  be  the  template  defined  on  the  m  x  n  array  X  by 

t(^i)={(^i).    (i  +  Hmod  m),j),    {i  -  l{mod  m),j),    {ij  +  l{mod  n)) 
(i,j  -  l{m.od  n))}. 

t  is  called  the  von  Neumann  averaging  (or  mean)  template.  This  template  can  also 
be  defined  by  the  rule 

r  1    ii\i-k\  +  \j-i\<i 

I,  0     otherwise 
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The  polynomial  associated  with  this  template,  via  the  isomorphism  6,  is 

p^{x,  2/)=  1  +  X  +  x"-l  +  ?/  +  y'"-l. 

Let  X  be  a  rectangular  array  with  m  rows  and  n  columns,  and  let  E  be  the 
template  defined  by  : 

E(x)  =  {(y,ex(y))  :  ex(x)  -  1   and     ex(y)  =  0     if  y  /  x}. 

E  is  the  identity  template  for  the  convolution  0;  that  is 

t®E  =  E®t  =  t. 

A  template  t  on  X  is  said  to  be  invertible  with  respect  to  ®  if  there  exists  a 
template  s  on  X  such  that 

tes  =  s©t  =  E. 


Proposition  6.4.  The  von  Neumann  averaging  template  t,  defined  on  a  finite 
rectangular  coordinate  set  with  m  rows  and  n  columns,  is  invertible  if  and  only  if 
Pl{ijjin:^n)  7^  0.  for  every  j  and  k. 

6.2.2.  The  Main  Theorem 

Theorem  6.5.  // X  is  an  m  x  m  rectangular  coordinate  set  in  Z  x  Z,  and  if  t 
denotes  the  von  Neumann  mean  filter  template  defined  on  X,  then  t  is  invertible  if 
and  only  if  neither  5  nor  6  divide  m. 
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ProoF: 

For  j,  A-  =  0, ...,  m  —  1  we  have 

Pt('^m,^m)  =  1  +  ^m  +  '^m    +  ^m  +  ^m 

=  1  +  2cos{2TTJlm)  +  2cos{2Trk/m). 
For  convenience,  we  will  write  a;  instead  of  Um- 
We  first  assume  that  one  of  the  integers  5  or  6  divides  m. 
Let  m  =  5J,  for  some  J.  If  we  let  j  =  d  and  k  =  2c?,  then 

p^  (a;-'  ,u-'j  =  l+a;+cj       +cj      +cj 


=  1  +  a;5  +  u?5     +  cjg  +  cj^ 

for  this  is  the  sum  of  the  five  fifth  roots  of  unity. 

Let  m  =  6c?,  for  some  d.  If  we  let  j  =  d  and  k  =  3d,  then 

=  1  +  2cos{Tr)  +  2cos{7r/3)  =  0. 

Hence  if  5  or  6  divides  m,  there  exist  two  integers  _;'  and  k  such  that  p^(u;5„,  w^) 
=  0  and,  by  Proposition  2.4,  the  template  t  is  not  invertible. 

WE  now  prove  the  sufficient  condition. 

Due  to  the  symmetry  oi  j  and  k,  in  the  polynomial  expression  p^(a;-^,u; '),  it  is 
enough  to  restrict  these  integers  to  0  <  j  <  A;  <  m/2. 

Let  j  be  the  minimal  value  such  that  p^(a;-',cj   )  =  0. 

Note  that  if 

q{x)  =  l  +  x^  +  x("'-^^J  +  x^-  +  .r('"-l)^', 
then  the  equation  p^(a7-^,a?  )  =  0,  is  equivalent  to  ^(u-')  =  0.  If  we  use  the  fact  that 
cyclotomic  polynomials  are  the  minimal  polynomials  of  the  primitive  roots  of  unity 
[3],  then  we  have  that  q(LO^)  =  0,  for  any  s  relatively  prime  to  m. 
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Claim:  li  p^{uj^ ,uj^)  =  0,  then  m/6  <  j  <  k  <  m/2. 

If  for  any  z  =  r*exp(m),  Arg{z)  =  a,  then  Arg{iu^)  =  Arg{u!''^).  Note  that 
in  order  to  find  two  integers  j  and  k  such  that  p^^{u)m,u>^)  =  0,  it  is  necessary 
that  7r/3  <  Arg{ijj^)  <  tt.  Since  Arg{u!^)  —  2iTJ/m,  the  previous  argument  gives 
'm/Q  ^  J  ^  m/2.  The  claim  is  proved. 

We  will  consider  two  cases:  j  divides  m  and  j  does  not  divide  m. 
Case  1:  j  divides  m. 

Using  the  claim,  we  know  that  2j  <  m  <  Qj.  If  j  divides  m,  then  m  is  of  the 
form  2j,  3i,  4j,  5j  or  6j. 

If  m  =  2j,  then  p^{lj^  .,lj  )  =  Oimplies  that  cos{2irk/m)  =  —1/2,  which  implies 

that  k  =  7n/3.  Since  k  is  an  integer,  3  divides  m.  Hence  6  divides  7n. 
If  m  =  3j,  then  p^(u;-^,i<;  )  =  0  implies  that  cos{2Trk/7n)  =  0,  which  implies 

that  k  =  m/6.  Since  3  and  4  divide  m,  6  divides  m. 
If  777  =  4j\  then  p^(a;-^,a;  )  =  0  implies  that  cos{2Trk/m,)  =  —1/2,  which  implies 

that  k  =  m/3.  Since  3  and  4  divide  m,  6  divides  tt?. 
If  m  =  5j  or  6j,  then  since  j  is  an  integer,  we  conclude  that  5  or  6  divides  m. 
Case  2:  j  does  not  divide  m.. 

Using  the  claim  again,  we  can  write  m  as  2j  +  r,  3j  +  r,  Aj  +  r,  or  5j  +  r  for 
some  nonzero  r  <  j. 

If  777  =  3j  +  r,  then  either  2  divides  ttj  or  2  does  not  divide  tt?..  In  the  case 
where  2  divides  m,  if  3  divides  ttt,  we  are  done.  Otherwise,  the  greatest 
commun  divisor  of  m  and  3,  denoted  by  gcd(m,3),  is  equal  to  1.  Thus 

q{u;^)  =  1  +  uj^i  +  uj-^J  +  u^^  +  u'^^  =  0. 
But  3j  —  m  —  r,  thus. 
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ui 


q{J)  =  1  +  u;'"-'-  +  a;-('"-'-^  +  u^^  + 

which  contradicts  the  minimality  of  j.  In  the  case  where  2  does  not  divide 
m,  then  gcd(?7i,4)  =  l,  and 

qiu"^)  =  1  +  u;4j  +  u;-4i  +  u^^  +  u'^^  =  0. 
But,  4 j  =  771  -  r  +  j,  thus, 

9(a;4)  =  1  +  u'^-r+J  +  ^-(m-r+j)  ^  ^4fc  ^  ^-4fc 

=  1  +  J-^  +  u;-(-^-'-^  +  u;4^'  +  u;-^^ 
which  contradicts  the  minimality  of  j. 
If  m  =  Aj  +  r,  then  either  3  divides  rri  or  3  does  not  divide  ttz.  In  the  case  where 
3  divides  772,  if  2  divides  777  we  are  done.  Otherwise,  gcd(77i,4)  =  l.  Thus, 
q{uj^)  =  1  +  L^4J  +  a;-4i  +  a;^^-  +  0;"^^  =  q. 

But,  4 j  =  777  —  r,  thus, 

/    4\        1    I       771-7-    ,   ,   -(777,-r)    I   ,  ,4fc    I   ,    -4A; 
q{ij  )  =  l+u;  +a;^         '+a;      +0; 

11       r    ,       — r    I      .5/;    ,    ,  ,—bk 

which  contradicts  the  minimality  of  j.  In  the  case  where  3  does  not  divide 
777,  then  gcd(r77.,9)  =  l  and 

q{u^)  =  1  +  a.9i  +  u-QJ  +  u^k  ^  ^-9k  ^  q 

But  9j  =  2777  -  2r  +  j,  thus, 

q{J)  =  1  +  u;2— 2'-+J  +  ^-(2m-2r+i)  ^  ^9^^  ^  ^-9^^ 
^  1  +  ^i-2r  +  ^-0--27-)  +  ^9^  +  ^-9^; 

which  contradicts  the  minimality  of  j,  since  |j  —  2r|  <  j. 
If  777  =  5j  +  r,  then  if  5  divides  777  we  are  done.  Otherwise,  gcd(77!,5)  =  1  and 

qiu^)  =  1  +  u^^  +  uJ-^^  +  oJ^''  +  ^-^^  =  0, 
since  uf*  is  a  primitive  777^^*  root  of  unity.  But  5j  =  m  -  r,  thus. 
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/     5\         1     I       m—r    ,       —(m—r)    ,       5fc    ,       —bk 
q[uj  )  =  l+u;  +  u)    ^         '  +  u      +  u 

It         V     I         — T     I         Ok     I         — ofc 
+  U      +  U)  +  iJ         +  tJ  , 

which  contradicts  the  minimality  of  j,  since  r  <  j.  This  concludes 
the  proof  of  Theorem  3.1 

Q.E.D. 


A  generalization  of  the  von  Neumann  averaging  template  is  the  one  defined  by 

a     if  r  =  I,  and  s  —  j 


t(fj)(^-5)   = 


b     otherwise 


t   = 


<a> 


where  a,  and  b  are  two  real  numbers. 

The  associated  polynomial  is  p{x,  y)  =  a  +  bx  +  bx^~    +  by  +  by^~  ,  and 

p{uj^,uj^)  =  a  +  bujj  +  buj~^  -\- bu^  +  bui~^ 
=  a  +  b{u^  +  uj~^  +  cj^  +  u~^). 

Proposition  6.6.    If  a  =  0,  then  the  template  t  is  invertible  if  and  only  if    m  is 
an  odd  number. 

Proof: 

Of  course,  we  will  assume,  here,  that  6  7^  0,  in  which  case  it  can  be  factored. 


If 


p{u^  ,U^)   =10^   +UJ     ^   +UJ^   +UJ      ^ 

=  2cos{27rj/m)  +  2cos{2Trk/m)  =  0, 
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then  the  equation  is  equivalent  to  cos{2Trj/m)    =    -cos{2nk/m),  which  implies 
2Trj/m  +  2Trk/m   =  tt,  and  therefore,  j  +  k  =  m/2.  Since  j  and  k  are  integers, 
this  proves  that  for  any  j,  and  k  =  m/2-j,  p{u\u^)  =  0.  Therefore,  the  template 
t  is  not  invertible  if  m  is  an  even  number. 

Q.E.D. 

Proposition  6.7.    If  a  ^  0  and  \a/b\  >  4,  then  the  template  t  is  invertible. 

Proof: 

If  pCu^^u;^)  =  g  +  bio^  +  6a;"-''  +  bu''  +  bu^"^ 

=  a  +  2b{cos{2Trj/m)  +  2cos{2Trk/m))  =  0, 

then  \cos{2Trj/m)  +  cos{2rk/m)\  =  \a/2b\  >  2,  which  is  a  contradiction.  Therefore, 
p{u^  ,uj^)  is  always  nonzero,  and  the  template  t  is  invertible. 

Q.E.D. 

6.3.      The  Characteristic  Polvnomial  Method 

Let  X  be  a  coordinate  set  with  m  rows  and  n  columns,  and  let  t  be  a  template 
defined  on  X.  Using  the  relationship  between  the  matrix  algebra  and  the  image 
algebra,  under  the  isomorphism  V',  we  know  that  the  corresponding  matrix  V'(t)  = 
Ml  is  an  mn  x  mn  block  matrix.  There  are  w?  blocks  and  each  block  is  an  n  x  7^ 
matrix. 

The  determinant  det(Mt  -  XImn),  where  Imn  is  the  unit  mn  x  mn  unit  or 
identity  matrix,  is  called  the  characteristic  polynomial  of  M^.  It  has  the  form 

p{\)  =  A"'"  +  a^„_iA"^"-l  +  ...  +  GiA  +  ao. 
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In  p{X),  the  constant  uq  is  the  determinant  of  M^,  while  the  coefficient  Grrm-l 
is  the  trace  of  M^  [12]. 

Theorem  6.8.   If  ao  is  nonzero,  then  M^^  is  invertible  and  its  inverse  is  given  by 
the  formula: 


M^-l  =  -i(M^"^"-^  +  arnn-lM^''~^  +  ...  +  a^In^n). 


Proof: 

By  the  Cayley-Hamilton  theorem,  we  know  that  the  matrix  M^  satisfies  its 
characteristic  polynomial  p(A),  i.e  p(Mt)  =  0.  Therefore, 

M(""  +  a^„_iM^""-^  +  ...  +  aiMt  +  ao/mn  =  0. 

The  previous  equation  can  be  rewritten  as 

M^"  +  a„,„_iM^"~^  +  ...  +  aiMt  =  -cto^^nn,  which  implies 

Mt[^(Mj""-^  +  a,„„_lMj""-2  +  ...  +  Gi7,nn)]  =  Imn- 

This  last  equation  gives 

Q.E.D. 

Let  E  be  the  identity  template  defined  in  Section  6.2.1,  and  let  t"  be  defined 

as  t  ©t  ©  ...®  t  n  times. 
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From  the  previous  theorem,  we  can  conclude  that,  given  a  shift  invariant  tem- 
plate t,  we  can  write  its  inverse  as 

t_i  =  zl(t-"-l  +  a,„„_2t"^"-2  +  ...  +  aiE). 
CO 

Since  the  inverse  of  a  block  toeplitz  matrix  with  toeplitz  blocks  is  not  neces- 
sarily a  matrix  with  the  same  properties,  we  see  that  even  the  inverse  of  a  very 
small  size  shift  invariant  operator  can  be  (and  most  probably  will  be)  a  full  size 
variant  operator.  The  previous  theorem  is  one  way  to  write  the  inverse  of  a  shift 
invariant  template  in  terms  of  the  template  itself  and  therefore,  in  terms  of  invariant 
templates. 

6.4.    Inverse  of  General  Shift  Invariant  Operators 

In  this  section,  we  discuss  some  properties  of  block  Toeplitz  matrices  with 
Toeplitz  blocks  as  they  are  directly  related  to  shift  invariant  templates  which  are, 
as  mentioned  previously,  used  to  implement  convolutions.  The  idea  here  is  to  apply 
some  results  on  Toeplitz  matrices  to  the  image  algebra. 

As  in  the  beginning  of  this  chapter,  the  matrices  of  this  section  are  indexed  by 
sets  of  the  form  {0,1,. ..,m  -  1}. 

Definition  6.9.  Let  M  =  {m.ij)  be  an  m  x  m.  matrix.  We  say  that  M  is  a 
Toeplitz  matrix  if  and  only  if  for  every  ij  G  {0, 1, ...,  m  -  1} ,  and  k  in  Z  such  that 
i  +  kj  +  k  e  {0, 1,  ...,m  -  1},  we  have  aij  =  aij^kj+k- 


Thus,  M  is  a  matrix  of  the  form 


a_l 


ao 


l-«-(m-l)         «-(7T!-2) 


..     a?n-2 
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A  Toeplitz  matrix  is  constant  along  any  diagonal  and  has  its  ij^"  entry  as  a 
function  of  {i  —  j)  rather  than  of  i  and  j  separately. 

Definition  6.10.  Let  M  —  (Mij)  be  an  m  x  m  block  matrix,  where  each  Mij  is 
an  n  X  n  matrix.  We  say  that  M  is  block  Toeplitz  with  Toeplitz  blocks  if  and  only  if 
M  is  of  the  form 


A-1 


Ao 


A-(m-l)      ^-(7n-2) 


Am-l 
Am-2 


where  for  every  i  in  {  —  {m  —  1),  (m  —  1)},  each  Aj  is  a  Toeplitz  matrix. 

The  following  theorem  was  proved  in  Gader  [8]. 

Theorem  6.11.  Let  X.  be  a  rectangular  coordinate  set.  If  t  is  a  template  on  X, 
and  M^  is  the  corresponding  matrix  under  ip,  then  M^  is  a  block  Toeplitz  matrix 
with  Toeplitz  blocks. 


Since  the  inverse  of  a  Toeplitz  matrix  is  not  necessarily  a  matrix  of  the  same 
form,  we  see  that  this  very  special  structure  of  Toeplitz  and  block  Toeplitz  matrices 
strongly  suggests  that  an  inversion  scheme  exploiting  this  structure  would  yield 
significant  savings  in  time  and  effort.  Several  researchers  have  tried  to  investigate 
Toeplitz  matrices  and  their  inverses.  Some  algorithms  like  the  Trench  algorithm 
proved  to  be  powerful  [33].    Other  methods  yielded  very  elegant  formulas  for  the 
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inverse  [1]. 


Let  us  first  give  some  facts  pertinent  only  to  Toeplitz  matrices. 

Let  M  be  an  m  X  m  Toeplitz  matrix;  u,  u,  two  m  x  1  vectors;  and  e*  the  transpose 


of  the  vector  (0,  ...0, 1,0,  ...0),  where  1  is  at  the  i^''  location. 


Fact  1: 


then  {um-l,Um-2.-,uo)M  =  (^m-l,  ^^771-2,  ■••,  ^^o)- 


Fact  2: 


If  Mu  =  eP,  and  Mv^e^, 
then  Um-i-q  -  Vm-\--p- 


Note  that  we  have 


\{Vm-\.-.nW\{u^,...,ih^-\f  =  (Vm-l,-,^0)[^^(^^0,-,^m-l)]- 


We  will  call  a  block  vector  of  dimension  m,  a  vector  in  which  all  entries  are 
n^n  matrices.  It  will  be  denoted  by  a  bold  letter.  We  will  denote  by  ej  the  vector 
(0, ...,  0,  /,  0, ...,  0),  where  the  n  x  n  unit  or  identity  matrix  I  is  at  the  i  entry,  and 
where  0  stands  for  the  n^n  zero  matrix.  1(^0,  ^1,  •••,  ^m-l)  represents  the  lower 
triangular  block  Toeplitz  matrix 


L  = 


^0 


0 


0 
0 


-1)      ^-(m-2)      •••      ^0- 
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In  a  manner  similar  to  the  lower  triangular  matrix,  we  denote  the  upper  trian- 
gular block  Toeplitz  matrix  by  U{Aq,  Ai, ...,  A^^i)  and  finally,  we  let  the  matrix  S 
denote  the  shift  matrix  S  =  L(0,  /,  0, ...,  0). 

Proposition  6.12.   Let  M  be  an  m  x  m  block  matrix  with  n  x  n  blocks.    If  there 
exist  column  block  vectors  x,  and  y  such  that 

Mx  =  eo   and  My  =  SMe^, 
then  M  is  invertible.  Moreover,  let  w  and  z  be  tivo  block  vectors  such  that 

wM  =  e(j    and  zMy  =  egMS, 
then 


M-1  =  AiBi- A2B2, 


where 

Ai  =  i:(xo,xi,...,x^_j  ,A2  =  i^(-yo,-yi,---,-ym-i)' 

Bi  =  t^(I, -zo,...,-z^_2)   ,B2  =  ?7(0,wo,...,w,„_2). 

In  the  implementation  of  this  theorem,  we  transpose  the  dreadful  problem  of 
inverting  an  mn  x  mn  block  Toeplitz  matrix  with  Toeplitz  blocks  to  the  problem  of 
solving  in  linear  systems  of  mn  equations  with  inn  variables. 


-73- 
In  the  application  of  this  theorem  to  image  algebra,  the  result  is  quite  inter- 
esting.  If  t  is  an  invariant  template  satisfying  the  conditions  of  the  theorem,  then 
there  exist  four  templates  ti,t2,t3,  and  t^  such  that 

t~^  =  ti  ®t2+t3et4, 

where  t^,  and  t^  are  shift  invariant  along  each  row,  and  t2  and  14  are  shift  invariant 
along  each  column 


CHAPTER  7 
CONCLUSION  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


This  dissertation  has  addressed  the  problems  of  template  decomposition  and 
template  inversion.  The  framework  is  based  on  the  image  algebra  and  its  relationship 
to  polynomial  and  matrix  algebra. 

In  particular,  it  was  shown  that: 

1)  A  two-dimensional  shift  invariant  operators  defined  on  a  rectangular  array 
can  be  decomposed  into  sums  and  products  of  smaller  size  operators.  The 
decomposition  methods  presented  here,  are  suitable  for  a  machine  with  a 
pipeline  architecture. 

2)  Two-dimensional  shift  invariant  operators  exhibiting  some  kind  of  symme- 
try require  a  much  smaller  number  of  operators  in  the  decomposition.  We 
focused  our  attention  on  this  type  of  operator  because  of  their  frequent  use 
in  image  processing. 

3)  The  techniques  were  extended  to  those  operators  which  may  not  be  shift 
invariant.  Necessary  and  sufficient  conditions  were  proved  to  test  for  the 
separability  of  a  variant  operator.  Necessary  and  sufficient  conditions  were 
proved  to  test  whether  or  not  a  variant  operator  can  be  written  as  the  sum 
of  separable  operators. 

4)  Necessary  and  sufficient  conditions  were  given  to  characterize  when  the  mean 
filter  based  on  the  von  Neumann  configuration  is  invertible.  Since  the  inverse 
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of  a  shift  invariant  operator  is  not  necessarily  shift  invariant  (actually,  there 
is  very  little  chance  it  is  shift  invariant),  methods  were  provided  for  writing 
the  inverse  in  ternns  of  shift  invariant  operators. 

We  now  list  some  suggestions  for  further  research. 

First,  generalize  the  results  obtained  in  Chapter  3  for  two-dimensional  shift 
invariant  operators  to  arbitrary  shift  invariant  operators  of  different  configuration. 
Circular  and  convex  configurations  would  be  of  particular  interest.  Can  the  poly- 
nomial decomposition  method  be  generalized  to  these  more  general  operators  ? 

The  decomposition  methods  for  variant  operators  would  certainely  be  imple- 
mented if  the  cost  in  computation  time  was  not  prohibitive.  Is  there  a  better  and 
more  direct  approach  to  their  decomposiiton  ?  Many  methods  given  in  Chapter 
4  involve  shift  invariant  and  variant  operators  in  the  decomposition  of  variant  op- 
erators. A  method  involving  only  shift  invariant  operators  will  surely  be  a  more 
efficient  one.  An  important  example  of  a  variant  template  is  the  Gabor  Transform. 
Can  the  results  presented  here  be  used  to  decompose  this  operator. 

A  generalization  of  the  condition  of  the  invertibility  of  the  von  Neumann  mean 
filter  defined  on  a  rectangular  array,  rather  than  on  a  square  array  should  be  proved 
in  the  the  very  near  future.  The  inverse  of  a  shift  invariant  operator  can  be  a  full  size 
variant  operator.  Methods  to  write  this  inverse  in  terms  of  an  acceptable  number 
of  shift  invariant  operators  can  find  their  way  in  areas  such  as  tomography,  remote 
sensing,  and  cryptography  among  others. 
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